简介本资源是一套面向光学工程、信息光学及计算成像方向初学者与高校实验教学的数字全息仿真实验MATLAB实现方案聚焦数字全息图生成、零级像抑制、波前再现等核心原理的编程验证。压缩包共2个文件1个BMP格式原始全息图数据、1个holographic.m主程序脚本总大小3.42MB结构精炼便于快速运行与代码剖析其中MATLAB脚本完整涵盖图像读取、傅里叶变换、空域高斯滤波去零级、逆变换及衍射再现全流程可直接用于课堂演示或课后复现。已有2423人学习下载适用于《信息光学》《计算全息》课程实验环节帮助学习者打通“光学原理—数值建模—MATLAB实现—图像分析”的完整链路切实掌握从干涉记录到三维物场重建的关键技术细节与调试逻辑。1. 从“拍立得”到“数字暗房”全息术的现代转身如果你对全息术的印象还停留在科幻电影里悬浮的立体影像或者博物馆里那些需要特定角度才能看到的静态图案那么是时候刷新一下认知了。传统的全息术就像一台复杂的“光学拍立得”它需要一束极其稳定的激光、一套精密的光学平台、一块特殊的感光干板以及一个暗房。整个过程对环境振动、空气流动都极为敏感制作成本高、周期长而且一旦拍坏就得从头再来。这极大地限制了全息术在科研和工程领域的普及与应用。而“数字全息仿真实验”则相当于把整个暗房搬进了电脑。我们不再需要真实的激光器和光学元件而是用MATLAB这样的数学工具在数字世界里模拟光的传播、干涉和衍射过程。这就像摄影师从胶片时代进入了数码时代拥有了一个功能强大的“数字暗房”。在这个暗房里你可以自由地调整“光源”的波长、角度改变“物体”的形状、位置甚至模拟现实中难以实现的理想条件而这一切只需要几行代码。对于光学工程、生物医学成像、微纳测量、信息安全等领域的研究者和学生来说这无疑打开了一扇低成本、高效率、可重复性极强的探索之门。本文的核心就是带你深入这个“数字暗房”手把手拆解如何用MATLAB实现一套完整的数字全息仿真流程。我们将从最基础的光波数学模型开始一步步构建出记录编码和再现解码全息图的全过程。你会发现那些看似神秘的全息原理其数学本质清晰而优雅。更重要的是我将分享在实际编码和调试过程中积累的诸多“坑点”和经验比如如何避免频谱混叠导致的图像失真如何选择最优的衍射算法来平衡精度与速度以及如何从再现的复杂光场中干净地分离出我们想要的物体像。无论你是刚接触计算光学的学生还是希望将仿真作为预研工具的工程师这篇超过五千字的深度指南都将为你提供一条从理论到实践的清晰路径。2. 数字全息的数学基石光场如何被“数字化”在动手写代码之前我们必须先理解数字全息仿真的核心——如何用离散的数学公式来描述连续的物理光场。这是所有后续工作的基础理解不透彻代码写出来也只是一知半解。2.1 光波的复数表达与采样定理在物理光学中一束单色光波在空间某一点 $(x, y, z)$ 和时间 $t$ 的振动可以用一个复数场 $U(x, y, z, t)$ 来表示。对于仿真我们通常关心的是空间分布且假设光波是单色的单一频率因此时间因子 $e^{-i\omega t}$ 可以分离出去我们只处理复振幅 $U(x, y, z)$。这个复振幅包含了光的振幅强度信息和相位波前形状信息这正是全息术能记录物体三维信息的根本。当我们用数字相机或计算机来模拟时面临的首要问题是将连续的空间 $(x, y)$ 离散化。假设我们在一个大小为 $L_x \times L_y$ 的平面上进行采样将其划分为 $N_x \times N_y$ 个像素点。那么采样间隔即像素尺寸为 $\Delta x L_x / N_x$ 和 $\Delta y L_y / N_y$。根据奈奎斯特采样定理要无失真地记录一个光场采样频率必须大于光场中最高空间频率的两倍。在数字全息中这直接决定了我们能模拟的物体的最大细节最高空间频率以及模拟区域的尺寸。如果采样不足就会发生频谱混叠在再现的图像中引入无法消除的噪声和伪影。注意在仿真中我们通常先确定采样点数 $N_x, N_y$ 和感兴趣的光学参数如波长 $\lambda$、传播距离 $d$然后反推所需的模拟区域尺寸 $L_x, L_y$以确保在整个计算过程中满足采样定理。这是一个关键的仿真设计步骤。2.2 干涉原理全息图的生成逻辑全息图的记录本质是干涉。一束光物光 $O$照射物体后其波前被物体调制携带了物体的信息。另一束光参考光 $R$作为基准直接照射记录介质。两者在记录平面相遇并发生干涉。干涉后的光强分布 $I(x, y)$ 为 $$I |O R|^2 |O|^2 |R|^2 O^R OR^$$ 其中$O^R$ 和 $OR^$ 是交叉项它们同时编码了物光的振幅和相位信息。$|O|^2$ 和 $|R|^2$ 是直流项通常构成背景噪声。数字全息仿真就是要在计算机中计算出这个 $I(x, y)$它就是我们模拟生成的“数字全息图”。在MATLAB中物光 $O$ 和参考光 $R$ 都是复数矩阵。例如一个常见的平面参考光可以表示为R A_R * exp(1i * 2*pi / lambda * (sin(theta_x)*X sin(theta_y)*Y));其中A_R是振幅theta_x和theta_y是参考光的倾斜角X和Y是坐标网格矩阵。物光 $O$ 的计算则更为复杂它需要模拟光从物体传播到记录平面的过程这就要用到下一节的核心——衍射计算。3. 光场传播模拟三种衍射算法的选择与实战如何计算光从物体平面传播到全息图平面或反过来的分布这是数字全息仿真的计算核心主要依靠标量衍射理论。最常用的有三种算法各有优劣选对算法直接影响仿真结果的准确性和计算效率。3.1 菲涅尔衍射Fresnel Diffraction与卷积法菲涅尔衍射积分是处理近场衍射的经典公式。其卷积形式为 $$U_2(x_2, y_2) \frac{e^{ikd}}{i\lambda d} \iint U_1(x_1, y_1) \exp\left( \frac{ik}{2d}[(x_2-x_1)^2(y_2-y_1)^2] \right) dx_1 dy_1$$ 可以看到输出场 $U_2$ 是输入场 $U_1$ 与一个二次相位因子菲涅尔传播子的卷积。在MATLAB中我们可以利用快速傅里叶变换FFT来实现卷积运算这就是卷积法Convolution Method。function U2 fresnel_convolution(U1, lambda, d, dx1, dy1) % U1: 输入平面复振幅 % lambda: 波长 % d: 传播距离 % dx1, dy1: 输入平面采样间隔 [M, N] size(U1); k 2*pi / lambda; % 生成输入平面的坐标 x1 (-N/2 : N/2-1) * dx1; y1 (-M/2 : M/2-1) * dy1; [X1, Y1] meshgrid(x1, y1); % 生成卷积核菲涅尔传播子 % 注意这里直接计算空间域的核对于大矩阵会非常慢且耗内存 % 因此通常采用基于FFT的频域相乘方式实现 % 以下为概念性代码展示原理 h exp(1i*k/(2*d) * (X1.^2 Y1.^2)); % 传播子 H fft2(ifftshift(h)); % 传播子的频域形式 U1_fft fft2(ifftshift(U1)); U2_fft U1_fft .* H; U2 fftshift(ifft2(U2_fft)); U2 U2 * exp(1i*k*d) / (1i*lambda*d); % 添加常数因子 end实操心得直接按上述公式实现卷积法需要计算一个非常大的卷积核效率低下。实际上更高效的方式是直接利用傅里叶变换的性质将空间域的卷积转换为频域的乘积。但即便如此标准的卷积法也存在一个关键限制输入和输出平面的采样间隔像素尺寸是相同的。这意味着如果你希望模拟的光传播距离d发生变化你无法同时自由控制输入和输出平面的物理尺寸这在实际仿真中很不方便。3.2 角谱法Angular Spectrum Method角谱法在频域处理衍射物理概念非常清晰。它将光场分解为不同方向传播的平面波角谱每个平面波在传播一段距离后仅产生一个相位延迟最后再合成新的光场。其数学表达为 $$U_2(x_2, y_2) \mathcal{F}^{-1} \left{ \mathcal{F}{U_1(x_1, y_1)} \cdot \exp\left( i k d \sqrt{1 - (\lambda f_x)^2 - (\lambda f_y)^2} \right) \right}$$ 其中 $\mathcal{F}$ 表示傅里叶变换$f_x, f_y$ 是空间频率坐标。角谱法的最大优点是精确。只要采样满足奈奎斯特条件它能精确满足亥姆霍兹方程理论上适用于任何传播距离从近场到远场。另一个巨大优势是输入和输出平面的采样间隔保持不变这给仿真带来了极大的便利。function U2 angular_spectrum(U1, lambda, d, dx, dy) % U1: 输入平面复振幅 % lambda: 波长 % d: 传播距离 % dx, dy: 采样间隔输入输出相同 [M, N] size(U1); k 2*pi / lambda; % 生成频域坐标 fx (-N/2 : N/2-1) / (N * dx); % 空间频率 (1/m) fy (-M/2 : M/2-1) / (M * dy); [FX, FY] meshgrid(fx, fy); % 计算传递函数 H exp(1i * k * d * sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); % 处理倏逝波空间频率过高根号内为负数的波 H(real(sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)) 0) 0; % 角谱传播 U1_fft fft2(ifftshift(U1)); U2_fft U1_fft .* ifftshift(H); % 注意传递函数需要ifftshift对齐 U2 fftshift(ifft2(U2_fft)); end踩坑记录角谱法中的传递函数H包含一个根号项。当 $(\lambda f_x)^2 (\lambda f_y)^2 1$ 时根号内为负数对应着倏逝波指数衰减不传播。在代码中必须处理这种情况通常将这部分传递函数设为0否则会引入数值不稳定。此外进行FFT和IFFT时fftshift和ifftshift的配对使用容易出错务必清楚每一步变换后数据的频率原点在矩阵的哪个位置。3.3 菲涅尔变换Fresnel Transform法当传播距离d满足菲涅尔近似条件时菲涅尔衍射积分可以简化为一个傅里叶变换的形式这就是菲涅尔变换法 $$U_2(x_2, y_2) \frac{e^{ikd}}{i\lambda d} e^{i\frac{k}{2d}(x_2^2y_2^2)} \cdot \mathcal{F} \left{ U_1(x_1, y_1) e^{i\frac{k}{2d}(x_1^2y_1^2)} \right}_{f_x\frac{x_2}{\lambda d}, f_y\frac{y_2}{\lambda d}}$$这个公式非常强大。它意味着输出平面的光场是输入平面光场乘上一个二次相位因子后的傅里叶变换再乘上另一个二次相位因子。更重要的是它揭示了输入和输出平面采样间隔的内在关系 $$\Delta x_2 \frac{\lambda d}{L_{x1}}, \quad \Delta y_2 \frac{\lambda d}{L_{y1}}$$ 其中 $L_{x1}, L_{y1}$ 是输入平面的物理尺寸。这意味着输出平面的像素尺寸 $\Delta x_2$ 由波长、传播距离和输入平面尺寸共同决定不是任意的。function [U2, x2, y2] fresnel_transform(U1, lambda, d, dx1, dy1) % U1: 输入平面复振幅 % lambda: 波长 % d: 传播距离 % dx1, dy1: 输入平面采样间隔 [M, N] size(U1); k 2*pi / lambda; % 输入平面尺寸和坐标 Lx1 N * dx1; Ly1 M * dy1; x1 (-N/2 : N/2-1) * dx1; y1 (-M/2 : M/2-1) * dy1; [X1, Y1] meshgrid(x1, y1); % 计算输出平面采样间隔和坐标 dx2 lambda * d / Lx1; dy2 lambda * d / Ly1; x2 (-N/2 : N/2-1) * dx2; y2 (-M/2 : M/2-1) * dy2; [X2, Y2] meshgrid(x2, y2); % 菲涅尔变换 phase_input exp(1i * k/(2*d) * (X1.^2 Y1.^2)); U1_modified U1 .* phase_input; U2_fft fft2(ifftshift(U1_modified)); phase_output exp(1i * k/(2*d) * (X2.^2 Y2.^2)); U2 (exp(1i*k*d) / (1i*lambda*d)) * phase_output .* fftshift(U2_fft); end算法选择指南算法优点缺点适用场景角谱法最精确输入输出采样间隔相同适用于任意距离计算量相对较大需处理倏逝波高精度仿真尤其是近场或需要固定像素尺寸的情况菲涅尔变换法计算速度快一次FFT公式简洁受菲涅尔近似条件限制输出像素尺寸被绑定满足菲涅尔近似的远场仿真快速预览卷积法概念直接计算效率低采样间隔固定教学演示原理实际工程中较少使用在我的大部分仿真工作中角谱法是首选因为它提供了最好的灵活性和精度。菲涅尔变换法则在需要快速计算、且传播距离较远时非常有用。4. 构建数字全息仿真系统从物体到全息图掌握了光场传播的工具后我们就可以搭建一个完整的数字全息仿真系统了。这个过程模拟了实际离轴全息记录的光路。4.1 物体建模与物光波前计算首先我们需要一个数字“物体”。最简单的模型是一个二维的振幅型物体如一个字母图案或相位型物体模拟一个微观的透明样本如细胞。更复杂的可以是三维点云或面型数据。% 示例创建一个振幅型物体字母‘F’ [M, N] deal(512, 512); % 采样点数 obj zeros(M, N); obj(200:300, 200:220) 1; % 竖线 obj(200:220, 200:300) 1; % 上横线 obj(260:280, 200:300) 1; % 中横线 % 物体通常位于某个平面z0我们赋予它一个初始的复振幅分布。 % 对于纯振幅物体相位可以设为0或一个随机相位模拟粗糙表面。 U_object obj; % 振幅为obj相位为0 % 或者加入随机相位模拟漫反射物体 % U_object obj .* exp(1i * 2*pi * rand(M, N));接下来模拟物光从物体传播到记录平面全息图平面的过程。假设记录平面距离物体平面为d_object。lambda 632.8e-9; % 氦氖激光波长单位米 d_object 0.1; % 传播距离0.1米 dx_obj 10e-6; % 物体平面采样间隔10微米 dy_obj 10e-6; % 使用角谱法计算记录平面上的物光波前 U_object_record angular_spectrum(U_object, lambda, d_object, dx_obj, dy_obj);这里U_object_record就是一个复数矩阵包含了物体信息调制后的光波的振幅和相位。4.2 参考光设计与全息图生成参考光通常设计为平面波或球面波。为了在再现时能将孪生像和零级像分离开我们采用离轴光路即让参考光以一定角度斜入射。% 生成记录平面的坐标网格 Lx_record N * dx_obj; % 注意角谱法传播后记录平面尺寸与物体平面相同 Ly_record M * dy_obj; x_record (-N/2 : N/2-1) * dx_obj; y_record (-M/2 : M/2-1) * dy_obj; [X_record, Y_record] meshgrid(x_record, y_record); % 设计一个倾斜的平面波作为参考光 theta_x_ref 0.5 * pi / 180; % 参考光在x方向倾斜0.5度弧度 theta_y_ref 0; % y方向不倾斜 k 2*pi / lambda; % 参考光复振幅 U_reference exp(1i * k * (sin(theta_x_ref)*X_record sin(theta_y_ref)*Y_record)); % 通常参考光强度远强于物光这里假设振幅为1。实际可调节比例。现在模拟干涉过程生成全息图光强分布。% 干涉光强 I_hologram abs(U_object_record U_reference).^2; % 此时 I_hologram 是一个实数矩阵模拟了CCD记录到的强度图。 % 为了模拟实际CCD的有限动态范围和量化噪声可以加入一些操作 % 1. 归一化到[0, 1] I_hologram I_hologram / max(I_hologram(:)); % 2. 模拟8位量化 I_hologram_quantized im2uint8(I_hologram); % 3. 加入高斯噪声模拟探测器噪声 I_hologram_noisy imnoise(I_hologram_quantized, gaussian, 0, 0.01); I_hologram double(I_hologram_noisy) / 255;生成的I_hologram就是我们的数字全息图。它可以被保存为图像文件模拟实际实验采集到的数据。4.3 全息图再现从编码光强到重建物体再现过程是记录的逆过程。我们用一束光再现光通常与参考光相同照射全息图然后计算衍射场最后在像平面得到重建的物体像。% 第一步用再现光照明全息图 % 假设使用与参考光相同的再现光这是最常见的共轭再现方式 U_illumination U_reference; % 再现光 % 全息图是光强需要将其视为一个振幅透射率或反射率调制 % 在线性记录条件下透射光场正比于光强乘以照明光 U_hologram I_hologram .* U_illumination; % 第二步将全息图平面视为新的光源向后传播到像平面 % 像平面通常位于原物体平面附近d_reconstruct ≈ -d_object d_reconstruct -d_object; % 负号表示反向传播 U_reconstructed angular_spectrum(U_hologram, lambda, d_reconstruct, dx_obj, dy_obj); % 第三步分析再现光场 % 再现光场 U_reconstructed 包含我们想要的原始物光波前虚像、 % 其共轭波实像孪生像以及零级衍射光直流项。 amplitude_recon abs(U_reconstructed); phase_recon angle(U_reconstructed); % 显示重建的振幅和相位 figure; subplot(1,2,1); imagesc(amplitude_recon); axis image; colormap gray; title(重建振幅); subplot(1,2,2); imagesc(phase_recon); axis image; colormap jet; title(重建相位);如果你直接运行上面的代码很可能会发现重建的图像一团糟几个像混叠在一起。这是因为我们还没有进行关键的频谱滤波步骤来分离出我们想要的像。5. 像的分离与优化频谱滤波与相位解包裹直接再现的全息图会同时产生零级像、原始像虚像和共轭像实像孪生像。为了得到清晰的物体像我们必须将它们分开。5.1 频域滤波在傅里叶空间“裁剪”目标离轴全息术的精妙之处在于通过引入参考光倾斜角使得物光信息被调制到了一个较高的载频上。在全息图的傅里叶频谱即角谱中零级和正负一级项对应原始像和共轭像在空间频率上是分离的。% 计算全息图的傅里叶频谱 I_fft fft2(ifftshift(I_hologram)); I_fft_shifted fftshift(I_fft); spectrum log(1 abs(I_fft_shifted)); % 对数显示增强对比 figure; imagesc(spectrum); colormap gray; axis image; title(全息图频谱); % 观察频谱你会看到中心亮斑零级和两侧对称的亮斑正负一级像我们的目标是通过一个滤波器只保留对应原始像的那一部分频谱。% 1. 设计一个带通滤波器以手动选择为例 [M, N] size(I_hologram); % 创建一个全零的滤波器矩阵 filter zeros(M, N); % 根据频谱观察确定原始像频谱中心的大致位置 (fx0, fy0) % 这需要根据参考光倾斜角计算或从频谱图中手动选取 fx0 50; % 示例在频域的x方向偏移 fy0 0; % y方向无偏移 radius 20; % 滤波器的半径 % 创建一个圆形掩模保留原始像频谱 [FX, FY] meshgrid(1:N, 1:M); filter(sqrt((FX - N/2 - fx0).^2 (FY - M/2 - fy0).^2) radius) 1; % 2. 应用滤波器 I_fft_filtered I_fft_shifted .* filter; % 3. 逆变换回空域得到滤波后的全息图 I_hologram_filtered fftshift(ifft2(ifftshift(I_fft_filtered))); % 4. 用滤波后的全息图进行再现 U_hologram_filtered I_hologram_filtered .* U_illumination; U_reconstructed_filtered angular_spectrum(U_hologram_filtered, lambda, d_reconstruct, dx_obj, dy_obj); amplitude_filtered abs(U_reconstructed_filtered);应用滤波后重建的振幅图像质量会显著提升孪生像和零级像的干扰被大幅抑制。经验技巧滤波器形状和大小是关键。圆形或矩形滤波器是常用选择。滤波器太小会损失物体高频细节图像模糊太大会引入其他级次的干扰。可以通过观察重建图像的质量来反复调整。更高级的方法是使用自动或半自动的频谱识别算法来定位滤波窗口。5.2 相位解包裹从缠绕相位到真实相位对于相位物体我们重建得到的是包裹相位Wrapped Phase其值被限制在 $[-\pi, \pi]$ 之间存在 $2\pi$ 的跳变。为了获得连续的、真实的物理相位分布需要进行相位解包裹Phase Unwrapping。% 重建的包裹相位 wrapped_phase angle(U_reconstructed_filtered); % 使用MATLAB内置的二维相位解包裹函数需要Image Processing Toolbox if license(test, image_toolbox) unwrapped_phase unwrap(wrapped_phase, [], 1); % 先解包裹行 unwrapped_phase unwrap(unwrapped_phase, [], 2); % 再解包裹列 else % 手动实现一个简单的质量引导路径积分法这里仅示意概念 % 实际应用推荐使用成熟的算法库如‘2D Phase Unwrapping’ by Miguel Arevallilo Herraez disp(Image Processing Toolbox未找到需手动实现或引入第三方解包裹算法。); end figure; subplot(1,2,1); imagesc(wrapped_phase); axis image; colormap jet; title(包裹相位); subplot(1,2,2); imagesc(unwrapped_phase); axis image; colormap jet; title(解包裹相位);避坑指南相位解包裹是数字全息定量测量中的难点。噪声、低调制区域振幅接近0或相位跳变过快超过采样定理都会导致解包裹错误产生“拉线”状的误差传播。在实际仿真中可以通过在物体前加入一个倾斜的相位平面相当于一个载频使相位变化更平缓从而降低解包裹难度。这类似于通信中的调制解调概念。6. 仿真进阶像差补偿与数值聚焦实际光学系统存在像差数字全息的优势在于可以在后处理中进行数字补偿。此外我们无需移动相机通过改变再现距离d_reconstruct就能实现数值聚焦这对三维成像至关重要。6.1 像差建模与数字补偿常见的像差如离焦、像散、彗差等都可以用泽尼克多项式来描述。我们可以在再现光中引入一个相反的相位板来抵消这些像差。% 假设我们已知系统存在主要的离焦像差和像散 % 在再现光场传播后引入一个补偿相位 [X_img, Y_img] meshgrid(x_record, y_record); % 像平面坐标 % 离焦系数 defocus_coeff 10; % 像散系数 astigmatism_coeff_x 5; astigmatism_coeff_y -5; % 构建像差相位 aberration_phase defocus_coeff * (X_img.^2 Y_img.^2) ... astigmatism_coeff_x * X_img.^2 astigmatism_coeff_y * Y_img.^2; % 补偿从重建光场中减去该像差相位 U_reconstructed_corrected U_reconstructed_filtered .* exp(-1i * aberration_phase);通过优化像差系数可以使重建图像的清晰度或条纹对比度达到最佳这个过程可以自动化例如使用图像清晰度评价函数如梯度平方和作为优化目标。6.2 数值聚焦与景深扩展对于稍厚的物体单一次重建只能对一个平面清晰成像。数字全息允许我们以不同的d_reconstruct反复进行再现计算从而得到一系列不同聚焦距离的切片图像实现三维层析。focus_range -0.11:0.001:-0.09; % 在物体平面前后各10mm范围内扫描 image_stack zeros(M, N, length(focus_range)); sharpness zeros(1, length(focus_range)); for idx 1:length(focus_range) d_focus focus_range(idx); U_temp angular_spectrum(U_hologram_filtered, lambda, d_focus, dx_obj, dy_obj); image_stack(:,:,idx) abs(U_temp).^2; % 保存强度图 % 计算该切片图像的清晰度例如用拉普拉斯算子的方差 lap del2(image_stack(:,:,idx)); sharpness(idx) var(lap(:)); end % 找到最清晰的聚焦平面 [~, best_idx] max(sharpness); best_focus_distance focus_range(best_idx); best_focused_image image_stack(:,:,best_idx); figure; plot(focus_range, sharpness, -o); xlabel(再现距离 (m)); ylabel(清晰度指标); title(自动聚焦曲线);这种方法可以自动确定物体的最佳聚焦位置并生成整个三维体积的数据对于生物细胞观测、微结构检测等应用非常有用。7. 从仿真到实践参数选择与性能优化心得经过上面几个章节一个完整的数字全息仿真框架已经建立。但在实际编码和调试中参数的选择和细节处理决定了仿真的成败与效率。这里分享一些关键的实战经验。1. 采样与尺寸的“先有鸡还是先有蛋”问题仿真开始时你需要确定采样点数N、像素尺寸dx、模拟区域尺寸L、波长λ和传播距离d。它们相互制约。我的建议是固定采样点数根据你的计算机内存和FFT效率先确定一个合适的N如512, 1024, 2048。2的幂次是FFT的最爱。由物理需求决定L或dx如果你关心物体的绝对尺寸就先定LL N * dx。如果你关心系统的分辨率由像素尺寸限制就先定dx。用角谱法验证对于角谱法只要dx λ/2满足采样定理且模拟区域L足够包含光场即可。一个快速检查方法是计算后传播的光场能量是否主要集中在模拟区域内边缘是否出现因周期边界假设导致的混叠。2. 参考光角度的“黄金分割”离轴角度θ的选择至关重要。角度太小频谱中各级像会混叠无法分离角度太大物光信息对应的空间频率可能超过探测器的奈奎斯特频率导致欠采样。理想情况是使物光频谱刚好与零频谱分开且留有足够的滤波窗口裕度。一个经验公式是载频 $f_c \sinθ / λ$ 应大于物体带宽B的1.5倍即 $f_c 1.5B$。在仿真中你可以先给一个小角度观察频谱再逐步调整。3. 计算速度与精度的权衡矩阵化操作避免在循环中进行逐像素计算充分利用MATLAB的矩阵运算和内置函数fft2,meshgrid,.^,.*等。选择合适的算法对于单次、高精度仿真用角谱法。如果需要快速扫描大量不同的再现距离如数值聚焦菲涅尔变换法可能更快因为其核心是一次FFT且改变距离d时只需重新计算输出平面的坐标和外部相位因子而角谱法则需要为每个d重新计算传递函数并做两次FFT。使用GPU加速如果仿真数据量很大如2048x2048以上且需要迭代优化可以考虑使用gpuArray将数据放到GPU上计算FFT在GPU上会有显著加速。4. 噪声与量化效应的模拟真实的实验全息图充满噪声。在仿真中加入噪声模型高斯噪声、泊松噪声和量化效应8位、12位量化能让你设计的重建算法更具鲁棒性。在滤波和相位解包裹步骤前对全息图进行适当的预处理如平场校正、减去背景的仿真也极具价值。数字全息仿真不仅仅是一个“正向建模”工具它更是一个强大的“逆向设计”和“算法验证”平台。你可以在其中自由地添加各种像差、噪声、振动模型然后测试你的补偿算法是否有效。你也可以用它来教学直观地展示“为什么参考光要倾斜”、“频谱滤波究竟做了什么”。当你真正吃透了这套仿真流程再去面对实际的实验数据和光路调试时你会拥有一种“上帝视角”般的透彻感。本文还有配套的精品资源点击获取 
长河君 · 2026/9/4 6:55:34 · 阅读 1.2k MATLAB数字全息仿真：从光场建模到图像重建的完整指南
简介本资源是一套面向光学工程、信息光学及计算成像方向初学者与高校实验教学的数字全息仿真实验MATLAB实现方案聚焦数字全息图生成、零级像抑制、波前再现等核心原理的编程验证。压缩包共2个文件1个BMP格式原始全息图数据、1个holographic.m主程序脚本总大小3.42MB结构精炼便于快速运行与代码剖析其中MATLAB脚本完整涵盖图像读取、傅里叶变换、空域高斯滤波去零级、逆变换及衍射再现全流程可直接用于课堂演示或课后复现。已有2423人学习下载适用于《信息光学》《计算全息》课程实验环节帮助学习者打通“光学原理—数值建模—MATLAB实现—图像分析”的完整链路切实掌握从干涉记录到三维物场重建的关键技术细节与调试逻辑。1. 从“拍立得”到“数字暗房”全息术的现代转身如果你对全息术的印象还停留在科幻电影里悬浮的立体影像或者博物馆里那些需要特定角度才能看到的静态图案那么是时候刷新一下认知了。传统的全息术就像一台复杂的“光学拍立得”它需要一束极其稳定的激光、一套精密的光学平台、一块特殊的感光干板以及一个暗房。整个过程对环境振动、空气流动都极为敏感制作成本高、周期长而且一旦拍坏就得从头再来。这极大地限制了全息术在科研和工程领域的普及与应用。而“数字全息仿真实验”则相当于把整个暗房搬进了电脑。我们不再需要真实的激光器和光学元件而是用MATLAB这样的数学工具在数字世界里模拟光的传播、干涉和衍射过程。这就像摄影师从胶片时代进入了数码时代拥有了一个功能强大的“数字暗房”。在这个暗房里你可以自由地调整“光源”的波长、角度改变“物体”的形状、位置甚至模拟现实中难以实现的理想条件而这一切只需要几行代码。对于光学工程、生物医学成像、微纳测量、信息安全等领域的研究者和学生来说这无疑打开了一扇低成本、高效率、可重复性极强的探索之门。本文的核心就是带你深入这个“数字暗房”手把手拆解如何用MATLAB实现一套完整的数字全息仿真流程。我们将从最基础的光波数学模型开始一步步构建出记录编码和再现解码全息图的全过程。你会发现那些看似神秘的全息原理其数学本质清晰而优雅。更重要的是我将分享在实际编码和调试过程中积累的诸多“坑点”和经验比如如何避免频谱混叠导致的图像失真如何选择最优的衍射算法来平衡精度与速度以及如何从再现的复杂光场中干净地分离出我们想要的物体像。无论你是刚接触计算光学的学生还是希望将仿真作为预研工具的工程师这篇超过五千字的深度指南都将为你提供一条从理论到实践的清晰路径。2. 数字全息的数学基石光场如何被“数字化”在动手写代码之前我们必须先理解数字全息仿真的核心——如何用离散的数学公式来描述连续的物理光场。这是所有后续工作的基础理解不透彻代码写出来也只是一知半解。2.1 光波的复数表达与采样定理在物理光学中一束单色光波在空间某一点 $(x, y, z)$ 和时间 $t$ 的振动可以用一个复数场 $U(x, y, z, t)$ 来表示。对于仿真我们通常关心的是空间分布且假设光波是单色的单一频率因此时间因子 $e^{-i\omega t}$ 可以分离出去我们只处理复振幅 $U(x, y, z)$。这个复振幅包含了光的振幅强度信息和相位波前形状信息这正是全息术能记录物体三维信息的根本。当我们用数字相机或计算机来模拟时面临的首要问题是将连续的空间 $(x, y)$ 离散化。假设我们在一个大小为 $L_x \times L_y$ 的平面上进行采样将其划分为 $N_x \times N_y$ 个像素点。那么采样间隔即像素尺寸为 $\Delta x L_x / N_x$ 和 $\Delta y L_y / N_y$。根据奈奎斯特采样定理要无失真地记录一个光场采样频率必须大于光场中最高空间频率的两倍。在数字全息中这直接决定了我们能模拟的物体的最大细节最高空间频率以及模拟区域的尺寸。如果采样不足就会发生频谱混叠在再现的图像中引入无法消除的噪声和伪影。注意在仿真中我们通常先确定采样点数 $N_x, N_y$ 和感兴趣的光学参数如波长 $\lambda$、传播距离 $d$然后反推所需的模拟区域尺寸 $L_x, L_y$以确保在整个计算过程中满足采样定理。这是一个关键的仿真设计步骤。2.2 干涉原理全息图的生成逻辑全息图的记录本质是干涉。一束光物光 $O$照射物体后其波前被物体调制携带了物体的信息。另一束光参考光 $R$作为基准直接照射记录介质。两者在记录平面相遇并发生干涉。干涉后的光强分布 $I(x, y)$ 为 $$I |O R|^2 |O|^2 |R|^2 O^R OR^$$ 其中$O^R$ 和 $OR^$ 是交叉项它们同时编码了物光的振幅和相位信息。$|O|^2$ 和 $|R|^2$ 是直流项通常构成背景噪声。数字全息仿真就是要在计算机中计算出这个 $I(x, y)$它就是我们模拟生成的“数字全息图”。在MATLAB中物光 $O$ 和参考光 $R$ 都是复数矩阵。例如一个常见的平面参考光可以表示为R A_R * exp(1i * 2*pi / lambda * (sin(theta_x)*X sin(theta_y)*Y));其中A_R是振幅theta_x和theta_y是参考光的倾斜角X和Y是坐标网格矩阵。物光 $O$ 的计算则更为复杂它需要模拟光从物体传播到记录平面的过程这就要用到下一节的核心——衍射计算。3. 光场传播模拟三种衍射算法的选择与实战如何计算光从物体平面传播到全息图平面或反过来的分布这是数字全息仿真的计算核心主要依靠标量衍射理论。最常用的有三种算法各有优劣选对算法直接影响仿真结果的准确性和计算效率。3.1 菲涅尔衍射Fresnel Diffraction与卷积法菲涅尔衍射积分是处理近场衍射的经典公式。其卷积形式为 $$U_2(x_2, y_2) \frac{e^{ikd}}{i\lambda d} \iint U_1(x_1, y_1) \exp\left( \frac{ik}{2d}[(x_2-x_1)^2(y_2-y_1)^2] \right) dx_1 dy_1$$ 可以看到输出场 $U_2$ 是输入场 $U_1$ 与一个二次相位因子菲涅尔传播子的卷积。在MATLAB中我们可以利用快速傅里叶变换FFT来实现卷积运算这就是卷积法Convolution Method。function U2 fresnel_convolution(U1, lambda, d, dx1, dy1) % U1: 输入平面复振幅 % lambda: 波长 % d: 传播距离 % dx1, dy1: 输入平面采样间隔 [M, N] size(U1); k 2*pi / lambda; % 生成输入平面的坐标 x1 (-N/2 : N/2-1) * dx1; y1 (-M/2 : M/2-1) * dy1; [X1, Y1] meshgrid(x1, y1); % 生成卷积核菲涅尔传播子 % 注意这里直接计算空间域的核对于大矩阵会非常慢且耗内存 % 因此通常采用基于FFT的频域相乘方式实现 % 以下为概念性代码展示原理 h exp(1i*k/(2*d) * (X1.^2 Y1.^2)); % 传播子 H fft2(ifftshift(h)); % 传播子的频域形式 U1_fft fft2(ifftshift(U1)); U2_fft U1_fft .* H; U2 fftshift(ifft2(U2_fft)); U2 U2 * exp(1i*k*d) / (1i*lambda*d); % 添加常数因子 end实操心得直接按上述公式实现卷积法需要计算一个非常大的卷积核效率低下。实际上更高效的方式是直接利用傅里叶变换的性质将空间域的卷积转换为频域的乘积。但即便如此标准的卷积法也存在一个关键限制输入和输出平面的采样间隔像素尺寸是相同的。这意味着如果你希望模拟的光传播距离d发生变化你无法同时自由控制输入和输出平面的物理尺寸这在实际仿真中很不方便。3.2 角谱法Angular Spectrum Method角谱法在频域处理衍射物理概念非常清晰。它将光场分解为不同方向传播的平面波角谱每个平面波在传播一段距离后仅产生一个相位延迟最后再合成新的光场。其数学表达为 $$U_2(x_2, y_2) \mathcal{F}^{-1} \left{ \mathcal{F}{U_1(x_1, y_1)} \cdot \exp\left( i k d \sqrt{1 - (\lambda f_x)^2 - (\lambda f_y)^2} \right) \right}$$ 其中 $\mathcal{F}$ 表示傅里叶变换$f_x, f_y$ 是空间频率坐标。角谱法的最大优点是精确。只要采样满足奈奎斯特条件它能精确满足亥姆霍兹方程理论上适用于任何传播距离从近场到远场。另一个巨大优势是输入和输出平面的采样间隔保持不变这给仿真带来了极大的便利。function U2 angular_spectrum(U1, lambda, d, dx, dy) % U1: 输入平面复振幅 % lambda: 波长 % d: 传播距离 % dx, dy: 采样间隔输入输出相同 [M, N] size(U1); k 2*pi / lambda; % 生成频域坐标 fx (-N/2 : N/2-1) / (N * dx); % 空间频率 (1/m) fy (-M/2 : M/2-1) / (M * dy); [FX, FY] meshgrid(fx, fy); % 计算传递函数 H exp(1i * k * d * sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)); % 处理倏逝波空间频率过高根号内为负数的波 H(real(sqrt(1 - (lambda*FX).^2 - (lambda*FY).^2)) 0) 0; % 角谱传播 U1_fft fft2(ifftshift(U1)); U2_fft U1_fft .* ifftshift(H); % 注意传递函数需要ifftshift对齐 U2 fftshift(ifft2(U2_fft)); end踩坑记录角谱法中的传递函数H包含一个根号项。当 $(\lambda f_x)^2 (\lambda f_y)^2 1$ 时根号内为负数对应着倏逝波指数衰减不传播。在代码中必须处理这种情况通常将这部分传递函数设为0否则会引入数值不稳定。此外进行FFT和IFFT时fftshift和ifftshift的配对使用容易出错务必清楚每一步变换后数据的频率原点在矩阵的哪个位置。3.3 菲涅尔变换Fresnel Transform法当传播距离d满足菲涅尔近似条件时菲涅尔衍射积分可以简化为一个傅里叶变换的形式这就是菲涅尔变换法 $$U_2(x_2, y_2) \frac{e^{ikd}}{i\lambda d} e^{i\frac{k}{2d}(x_2^2y_2^2)} \cdot \mathcal{F} \left{ U_1(x_1, y_1) e^{i\frac{k}{2d}(x_1^2y_1^2)} \right}_{f_x\frac{x_2}{\lambda d}, f_y\frac{y_2}{\lambda d}}$$这个公式非常强大。它意味着输出平面的光场是输入平面光场乘上一个二次相位因子后的傅里叶变换再乘上另一个二次相位因子。更重要的是它揭示了输入和输出平面采样间隔的内在关系 $$\Delta x_2 \frac{\lambda d}{L_{x1}}, \quad \Delta y_2 \frac{\lambda d}{L_{y1}}$$ 其中 $L_{x1}, L_{y1}$ 是输入平面的物理尺寸。这意味着输出平面的像素尺寸 $\Delta x_2$ 由波长、传播距离和输入平面尺寸共同决定不是任意的。function [U2, x2, y2] fresnel_transform(U1, lambda, d, dx1, dy1) % U1: 输入平面复振幅 % lambda: 波长 % d: 传播距离 % dx1, dy1: 输入平面采样间隔 [M, N] size(U1); k 2*pi / lambda; % 输入平面尺寸和坐标 Lx1 N * dx1; Ly1 M * dy1; x1 (-N/2 : N/2-1) * dx1; y1 (-M/2 : M/2-1) * dy1; [X1, Y1] meshgrid(x1, y1); % 计算输出平面采样间隔和坐标 dx2 lambda * d / Lx1; dy2 lambda * d / Ly1; x2 (-N/2 : N/2-1) * dx2; y2 (-M/2 : M/2-1) * dy2; [X2, Y2] meshgrid(x2, y2); % 菲涅尔变换 phase_input exp(1i * k/(2*d) * (X1.^2 Y1.^2)); U1_modified U1 .* phase_input; U2_fft fft2(ifftshift(U1_modified)); phase_output exp(1i * k/(2*d) * (X2.^2 Y2.^2)); U2 (exp(1i*k*d) / (1i*lambda*d)) * phase_output .* fftshift(U2_fft); end算法选择指南算法优点缺点适用场景角谱法最精确输入输出采样间隔相同适用于任意距离计算量相对较大需处理倏逝波高精度仿真尤其是近场或需要固定像素尺寸的情况菲涅尔变换法计算速度快一次FFT公式简洁受菲涅尔近似条件限制输出像素尺寸被绑定满足菲涅尔近似的远场仿真快速预览卷积法概念直接计算效率低采样间隔固定教学演示原理实际工程中较少使用在我的大部分仿真工作中角谱法是首选因为它提供了最好的灵活性和精度。菲涅尔变换法则在需要快速计算、且传播距离较远时非常有用。4. 构建数字全息仿真系统从物体到全息图掌握了光场传播的工具后我们就可以搭建一个完整的数字全息仿真系统了。这个过程模拟了实际离轴全息记录的光路。4.1 物体建模与物光波前计算首先我们需要一个数字“物体”。最简单的模型是一个二维的振幅型物体如一个字母图案或相位型物体模拟一个微观的透明样本如细胞。更复杂的可以是三维点云或面型数据。% 示例创建一个振幅型物体字母‘F’ [M, N] deal(512, 512); % 采样点数 obj zeros(M, N); obj(200:300, 200:220) 1; % 竖线 obj(200:220, 200:300) 1; % 上横线 obj(260:280, 200:300) 1; % 中横线 % 物体通常位于某个平面z0我们赋予它一个初始的复振幅分布。 % 对于纯振幅物体相位可以设为0或一个随机相位模拟粗糙表面。 U_object obj; % 振幅为obj相位为0 % 或者加入随机相位模拟漫反射物体 % U_object obj .* exp(1i * 2*pi * rand(M, N));接下来模拟物光从物体传播到记录平面全息图平面的过程。假设记录平面距离物体平面为d_object。lambda 632.8e-9; % 氦氖激光波长单位米 d_object 0.1; % 传播距离0.1米 dx_obj 10e-6; % 物体平面采样间隔10微米 dy_obj 10e-6; % 使用角谱法计算记录平面上的物光波前 U_object_record angular_spectrum(U_object, lambda, d_object, dx_obj, dy_obj);这里U_object_record就是一个复数矩阵包含了物体信息调制后的光波的振幅和相位。4.2 参考光设计与全息图生成参考光通常设计为平面波或球面波。为了在再现时能将孪生像和零级像分离开我们采用离轴光路即让参考光以一定角度斜入射。% 生成记录平面的坐标网格 Lx_record N * dx_obj; % 注意角谱法传播后记录平面尺寸与物体平面相同 Ly_record M * dy_obj; x_record (-N/2 : N/2-1) * dx_obj; y_record (-M/2 : M/2-1) * dy_obj; [X_record, Y_record] meshgrid(x_record, y_record); % 设计一个倾斜的平面波作为参考光 theta_x_ref 0.5 * pi / 180; % 参考光在x方向倾斜0.5度弧度 theta_y_ref 0; % y方向不倾斜 k 2*pi / lambda; % 参考光复振幅 U_reference exp(1i * k * (sin(theta_x_ref)*X_record sin(theta_y_ref)*Y_record)); % 通常参考光强度远强于物光这里假设振幅为1。实际可调节比例。现在模拟干涉过程生成全息图光强分布。% 干涉光强 I_hologram abs(U_object_record U_reference).^2; % 此时 I_hologram 是一个实数矩阵模拟了CCD记录到的强度图。 % 为了模拟实际CCD的有限动态范围和量化噪声可以加入一些操作 % 1. 归一化到[0, 1] I_hologram I_hologram / max(I_hologram(:)); % 2. 模拟8位量化 I_hologram_quantized im2uint8(I_hologram); % 3. 加入高斯噪声模拟探测器噪声 I_hologram_noisy imnoise(I_hologram_quantized, gaussian, 0, 0.01); I_hologram double(I_hologram_noisy) / 255;生成的I_hologram就是我们的数字全息图。它可以被保存为图像文件模拟实际实验采集到的数据。4.3 全息图再现从编码光强到重建物体再现过程是记录的逆过程。我们用一束光再现光通常与参考光相同照射全息图然后计算衍射场最后在像平面得到重建的物体像。% 第一步用再现光照明全息图 % 假设使用与参考光相同的再现光这是最常见的共轭再现方式 U_illumination U_reference; % 再现光 % 全息图是光强需要将其视为一个振幅透射率或反射率调制 % 在线性记录条件下透射光场正比于光强乘以照明光 U_hologram I_hologram .* U_illumination; % 第二步将全息图平面视为新的光源向后传播到像平面 % 像平面通常位于原物体平面附近d_reconstruct ≈ -d_object d_reconstruct -d_object; % 负号表示反向传播 U_reconstructed angular_spectrum(U_hologram, lambda, d_reconstruct, dx_obj, dy_obj); % 第三步分析再现光场 % 再现光场 U_reconstructed 包含我们想要的原始物光波前虚像、 % 其共轭波实像孪生像以及零级衍射光直流项。 amplitude_recon abs(U_reconstructed); phase_recon angle(U_reconstructed); % 显示重建的振幅和相位 figure; subplot(1,2,1); imagesc(amplitude_recon); axis image; colormap gray; title(重建振幅); subplot(1,2,2); imagesc(phase_recon); axis image; colormap jet; title(重建相位);如果你直接运行上面的代码很可能会发现重建的图像一团糟几个像混叠在一起。这是因为我们还没有进行关键的频谱滤波步骤来分离出我们想要的像。5. 像的分离与优化频谱滤波与相位解包裹直接再现的全息图会同时产生零级像、原始像虚像和共轭像实像孪生像。为了得到清晰的物体像我们必须将它们分开。5.1 频域滤波在傅里叶空间“裁剪”目标离轴全息术的精妙之处在于通过引入参考光倾斜角使得物光信息被调制到了一个较高的载频上。在全息图的傅里叶频谱即角谱中零级和正负一级项对应原始像和共轭像在空间频率上是分离的。% 计算全息图的傅里叶频谱 I_fft fft2(ifftshift(I_hologram)); I_fft_shifted fftshift(I_fft); spectrum log(1 abs(I_fft_shifted)); % 对数显示增强对比 figure; imagesc(spectrum); colormap gray; axis image; title(全息图频谱); % 观察频谱你会看到中心亮斑零级和两侧对称的亮斑正负一级像我们的目标是通过一个滤波器只保留对应原始像的那一部分频谱。% 1. 设计一个带通滤波器以手动选择为例 [M, N] size(I_hologram); % 创建一个全零的滤波器矩阵 filter zeros(M, N); % 根据频谱观察确定原始像频谱中心的大致位置 (fx0, fy0) % 这需要根据参考光倾斜角计算或从频谱图中手动选取 fx0 50; % 示例在频域的x方向偏移 fy0 0; % y方向无偏移 radius 20; % 滤波器的半径 % 创建一个圆形掩模保留原始像频谱 [FX, FY] meshgrid(1:N, 1:M); filter(sqrt((FX - N/2 - fx0).^2 (FY - M/2 - fy0).^2) radius) 1; % 2. 应用滤波器 I_fft_filtered I_fft_shifted .* filter; % 3. 逆变换回空域得到滤波后的全息图 I_hologram_filtered fftshift(ifft2(ifftshift(I_fft_filtered))); % 4. 用滤波后的全息图进行再现 U_hologram_filtered I_hologram_filtered .* U_illumination; U_reconstructed_filtered angular_spectrum(U_hologram_filtered, lambda, d_reconstruct, dx_obj, dy_obj); amplitude_filtered abs(U_reconstructed_filtered);应用滤波后重建的振幅图像质量会显著提升孪生像和零级像的干扰被大幅抑制。经验技巧滤波器形状和大小是关键。圆形或矩形滤波器是常用选择。滤波器太小会损失物体高频细节图像模糊太大会引入其他级次的干扰。可以通过观察重建图像的质量来反复调整。更高级的方法是使用自动或半自动的频谱识别算法来定位滤波窗口。5.2 相位解包裹从缠绕相位到真实相位对于相位物体我们重建得到的是包裹相位Wrapped Phase其值被限制在 $[-\pi, \pi]$ 之间存在 $2\pi$ 的跳变。为了获得连续的、真实的物理相位分布需要进行相位解包裹Phase Unwrapping。% 重建的包裹相位 wrapped_phase angle(U_reconstructed_filtered); % 使用MATLAB内置的二维相位解包裹函数需要Image Processing Toolbox if license(test, image_toolbox) unwrapped_phase unwrap(wrapped_phase, [], 1); % 先解包裹行 unwrapped_phase unwrap(unwrapped_phase, [], 2); % 再解包裹列 else % 手动实现一个简单的质量引导路径积分法这里仅示意概念 % 实际应用推荐使用成熟的算法库如‘2D Phase Unwrapping’ by Miguel Arevallilo Herraez disp(Image Processing Toolbox未找到需手动实现或引入第三方解包裹算法。); end figure; subplot(1,2,1); imagesc(wrapped_phase); axis image; colormap jet; title(包裹相位); subplot(1,2,2); imagesc(unwrapped_phase); axis image; colormap jet; title(解包裹相位);避坑指南相位解包裹是数字全息定量测量中的难点。噪声、低调制区域振幅接近0或相位跳变过快超过采样定理都会导致解包裹错误产生“拉线”状的误差传播。在实际仿真中可以通过在物体前加入一个倾斜的相位平面相当于一个载频使相位变化更平缓从而降低解包裹难度。这类似于通信中的调制解调概念。6. 仿真进阶像差补偿与数值聚焦实际光学系统存在像差数字全息的优势在于可以在后处理中进行数字补偿。此外我们无需移动相机通过改变再现距离d_reconstruct就能实现数值聚焦这对三维成像至关重要。6.1 像差建模与数字补偿常见的像差如离焦、像散、彗差等都可以用泽尼克多项式来描述。我们可以在再现光中引入一个相反的相位板来抵消这些像差。% 假设我们已知系统存在主要的离焦像差和像散 % 在再现光场传播后引入一个补偿相位 [X_img, Y_img] meshgrid(x_record, y_record); % 像平面坐标 % 离焦系数 defocus_coeff 10; % 像散系数 astigmatism_coeff_x 5; astigmatism_coeff_y -5; % 构建像差相位 aberration_phase defocus_coeff * (X_img.^2 Y_img.^2) ... astigmatism_coeff_x * X_img.^2 astigmatism_coeff_y * Y_img.^2; % 补偿从重建光场中减去该像差相位 U_reconstructed_corrected U_reconstructed_filtered .* exp(-1i * aberration_phase);通过优化像差系数可以使重建图像的清晰度或条纹对比度达到最佳这个过程可以自动化例如使用图像清晰度评价函数如梯度平方和作为优化目标。6.2 数值聚焦与景深扩展对于稍厚的物体单一次重建只能对一个平面清晰成像。数字全息允许我们以不同的d_reconstruct反复进行再现计算从而得到一系列不同聚焦距离的切片图像实现三维层析。focus_range -0.11:0.001:-0.09; % 在物体平面前后各10mm范围内扫描 image_stack zeros(M, N, length(focus_range)); sharpness zeros(1, length(focus_range)); for idx 1:length(focus_range) d_focus focus_range(idx); U_temp angular_spectrum(U_hologram_filtered, lambda, d_focus, dx_obj, dy_obj); image_stack(:,:,idx) abs(U_temp).^2; % 保存强度图 % 计算该切片图像的清晰度例如用拉普拉斯算子的方差 lap del2(image_stack(:,:,idx)); sharpness(idx) var(lap(:)); end % 找到最清晰的聚焦平面 [~, best_idx] max(sharpness); best_focus_distance focus_range(best_idx); best_focused_image image_stack(:,:,best_idx); figure; plot(focus_range, sharpness, -o); xlabel(再现距离 (m)); ylabel(清晰度指标); title(自动聚焦曲线);这种方法可以自动确定物体的最佳聚焦位置并生成整个三维体积的数据对于生物细胞观测、微结构检测等应用非常有用。7. 从仿真到实践参数选择与性能优化心得经过上面几个章节一个完整的数字全息仿真框架已经建立。但在实际编码和调试中参数的选择和细节处理决定了仿真的成败与效率。这里分享一些关键的实战经验。1. 采样与尺寸的“先有鸡还是先有蛋”问题仿真开始时你需要确定采样点数N、像素尺寸dx、模拟区域尺寸L、波长λ和传播距离d。它们相互制约。我的建议是固定采样点数根据你的计算机内存和FFT效率先确定一个合适的N如512, 1024, 2048。2的幂次是FFT的最爱。由物理需求决定L或dx如果你关心物体的绝对尺寸就先定LL N * dx。如果你关心系统的分辨率由像素尺寸限制就先定dx。用角谱法验证对于角谱法只要dx λ/2满足采样定理且模拟区域L足够包含光场即可。一个快速检查方法是计算后传播的光场能量是否主要集中在模拟区域内边缘是否出现因周期边界假设导致的混叠。2. 参考光角度的“黄金分割”离轴角度θ的选择至关重要。角度太小频谱中各级像会混叠无法分离角度太大物光信息对应的空间频率可能超过探测器的奈奎斯特频率导致欠采样。理想情况是使物光频谱刚好与零频谱分开且留有足够的滤波窗口裕度。一个经验公式是载频 $f_c \sinθ / λ$ 应大于物体带宽B的1.5倍即 $f_c 1.5B$。在仿真中你可以先给一个小角度观察频谱再逐步调整。3. 计算速度与精度的权衡矩阵化操作避免在循环中进行逐像素计算充分利用MATLAB的矩阵运算和内置函数fft2,meshgrid,.^,.*等。选择合适的算法对于单次、高精度仿真用角谱法。如果需要快速扫描大量不同的再现距离如数值聚焦菲涅尔变换法可能更快因为其核心是一次FFT且改变距离d时只需重新计算输出平面的坐标和外部相位因子而角谱法则需要为每个d重新计算传递函数并做两次FFT。使用GPU加速如果仿真数据量很大如2048x2048以上且需要迭代优化可以考虑使用gpuArray将数据放到GPU上计算FFT在GPU上会有显著加速。4. 噪声与量化效应的模拟真实的实验全息图充满噪声。在仿真中加入噪声模型高斯噪声、泊松噪声和量化效应8位、12位量化能让你设计的重建算法更具鲁棒性。在滤波和相位解包裹步骤前对全息图进行适当的预处理如平场校正、减去背景的仿真也极具价值。数字全息仿真不仅仅是一个“正向建模”工具它更是一个强大的“逆向设计”和“算法验证”平台。你可以在其中自由地添加各种像差、噪声、振动模型然后测试你的补偿算法是否有效。你也可以用它来教学直观地展示“为什么参考光要倾斜”、“频谱滤波究竟做了什么”。当你真正吃透了这套仿真流程再去面对实际的实验数据和光路调试时你会拥有一种“上帝视角”般的透彻感。本文还有配套的精品资源点击获取