企业官网建设流程全解析

简介：这是一份基于相干衍射成像（CDI）的MATLAB模拟实现源码包，面向光学成像、计算成像以及信号处理方向的研究生、科研人员与工程师，尤其适合刚接触相干衍射成像、希望通过Matlab快速建立模拟链路的学习者。资源包内含2个m文件，均为MATLAB函数/脚本，从波前传播到迭代重建的关键步骤均有体现；压缩包整体仅1KB，轻量精简，便于直接阅读、调试与二次修改。目前已有223人学习/下载，对于入门相干衍射成像仿真的参考价值较为明确。通过该源码，可了解衍射传播模型的构建方式、迭代重建过程的循环策略与数据流组织思路；同时，源码结构清晰，适合在此基础上发展自定义相位恢复或成像优化算法，节省从零搭建仿真环境的时间，尤其适合作为课程设计、课题预研或论文复现的起步代码。

1. 只拍到强度却要重建复振幅：CDI 模拟到底在造什么数据

光学成像里最反直觉的一件事，是你丢掉了相位还能把物体还原出来。相干衍射成像（CDI）记录的是远场衍射强度，相位信息在探测器上根本没被存下来，但通过对物体施加“有限支撑”这个先验，再用交替投影迭代，居然能同时恢复振幅和相位。模拟源码的价值，就是把实验里最难控制的照明分布、噪声水平、支撑域误差全换成已知量，让你在 MATLAB 里一遍遍重跑同一条迭代链路，观察哪个参数真正影响收敛。这篇博文面向用 MATLAB 做光学仿真、数字全息、相衬成像或刚接触相位恢复的工程师，用一个最小 CDI 模拟拆解源码里的关键环节：物理模型怎么落成 fft2、支撑域怎么初始化、HIO 迭代为什么在某个参数下停住，以及最后怎么确认重建结果是可信的。

2. 相干衍射成像模拟的物理模型与 MATLAB 矩阵化

2.1 从夫琅禾费衍射到二维 FFT：模拟源码的第一行等式

CDI 的前向模型很干净：物体被相干光照明后，出射波在远场传播，探测器平面上的复振幅等于物面出射波的傅里叶变换。写成离散形式就是一句fft2，这也是几乎所有 MATLAB 模拟源码的第一行核心运算。

I(u,v) = | FFT{ P(x,y) · O(x,y) } |²

其中O(x,y)是物体的复振幅透过率，P(x,y)是照明光斑，整个乘积就是入射波与物体作用后的出射波。注意这里有一个新手容易忽略、老手也容易踩的点：波长、距离、像素物理尺寸这些参数在离散 FFT 模型里不会显式出现，因为它们只影响频域坐标的缩放比例，不影响重建算法本身。模拟里把波长距离归一化，矩阵索引差 1 格，就相当于物理空间里差一个固定尺度。

MATLAB 里要写对这一步，关键不是公式，而是象限原点。fft2默认原点在矩阵左上角，物理上我们希望波前原点在矩阵中心，所以标准写法是输入先做ifftshift，输出再做fftshift。这个配对写反一次，重建出来的物体就会带一个线性相位斜坡，看起来像整体偏转了。

wave = obj .* illum; % 出射波场 F = fftshift(fft2(ifftshift(wave))); % 中心化二维FFT I = abs(F).^2; % 探测器强度

这里ifftshift把矩阵中心移到(1,1)供fft2处理，fftshift再把频域结果恢复到以中心为原点的布局。后续做逆变换时要完全反向配对，否则空间域和频域的原点就对不齐，自相关支撑也会算歪。

2.2 支撑域与过采样比：为什么矩阵尺寸不能随便给

CDI 能成立，不是靠算法多聪明，而是靠一个物理先验：物体占据了有限大小的区域，称为支撑域（support）。支撑域外的空间里，出射波必须严格为零。交替投影算法做的事，就是让当前估计反复在两个约束集合之间来回投影：傅里叶域要求振幅等于测量的sqrt(I)，实空间要求支撑域外为零。两个约束的交集越来越小，解才越来越唯一。

这就引出一个参数：过采样比。它定义为探测器边长与物体支撑宽度的比值，一维物体理论下限是 2，二维实际经验也要取到 2 到 4。模拟里最常犯的错误是把物体铺满了整个矩阵，比如坐落在 256×256 矩阵里的一个 200×200 矩形块。这种配置下，缺失相位信息的自由度太多，相位恢复算法基本不会收敛，误差曲线会从第 10 次迭代开始就平行下滑，再跑 2000 轮也没用。

实用的生成规则是：物体直径控制在矩阵边长的 0.25 到 0.4 倍，让物体周围留出至少一倍直径的零背景。后面章节给出的模拟参数会沿用这个经验值。

提示：过采样比不是越高越好。支撑域只占矩阵的 1/16 时，频域采样确实更密，但有效信号能量占比太低，噪声影响会被放大。0.3 倍边长附近是模拟里比较稳的区间。

2.3 把探测器参数映射成 MATLAB 变量：一张参数表

实验里的探测器参数看起来和矩阵运算隔得很远，但模拟源码里每一个都能对应到一个变量或一行处理。下面这张表是我在模拟里最常用的映射关系，也方便把实验导出的数据套进同一套脚本。

实验里从相机导出的强度图通常是 CSV 或 TIFF 格式，CSV 场景下用readmatrix('data.csv')读进来就是矩阵，后面的流程和模拟完全一致。唯一要确认的是数值范围和单位，模拟代码里一般会先做一次归一化，保证强度峰值不是几万量级的大数，避免后续poissrnd或sqrt出现数值溢出。

3. 用 MATLAB 源码搭一套 CDI 模拟的最小可行流程

3.1 生成复振幅物体与照明光斑

模拟的第一步是制造一个“已知真值”的复振幅物体。振幅部分可以用图像处理工具箱里的phantom生成一个类似组织结构的灰度图，相位部分叠加一个平滑面形，让重建任务同时包含振幅对比和相位延迟。

% CDI 模拟参数区 N = 256; % 探测器像素数（矩阵边长） obj_r = round(N * 0.18); % 物体半径，过采样比约 2.8 rng(0); % 固定随机种子，保证可复现 % 生成复振幅物体：幅值 + 相位 amp = phantom(N); % Shepp-Logan 模型作幅值 [X, Y] = meshgrid(1:N, 1:N); phase = 0.4 * peaks(N); % 平滑相位面形，峰值约 2.4 rad obj = amp .* exp(1i * phase); % 复数透过率 % 照明光斑：平面波（全1）；高斯照明时取消注释下一段 illum = ones(N); % illum = exp(-((X-N/2).^2 + (Y-N/2).^2) / (2*(0.5*obj_r)^2)); wave = obj .* illum; % 出射波场

phantom的输出范围是 0 到 1，peaks函数输出大约在 -6 到 6 之间，乘 0.4 后相位幅度约 ±2.4 rad，既有明显的相位包裹，又不至于让相位梯度大得离谱。物体半径取0.18 * N，物体外圈有超过自身直径两倍的零背景，过采样比满足要求。

照明部分默认用平面波。想模拟聚焦照明时，把注释掉的高斯光斑换上去，高斯宽度用物体半径的一半，重建时支撑域会自动适应照明轮廓，这也是模拟相对实验的便捷之处。

3.2 前向传播与强度记录

前向传播只需要一次中心化 FFT，然后取模平方。为了接近真实实验，加泊松噪声模拟光子计数统计涨落，并加一个中心遮挡来模拟光束阻挡。

F = fftshift(fft2(ifftshift(wave))); % 远场衍射 I = abs(F).^2; % 泊松噪声：控制峰值光子数，影响信噪比 photon_count = 1e4; % 峰值光子数 I_norm = I / max(I(:)) * photon_count; I_noisy = poissrnd(I_norm); I_meas = I_noisy / photon_count * max(I(:)); % 还原到原始量级 % 中心光束阻挡（beamstop） r_bs = 3; [xx, yy] = meshgrid(1:N, 1:N); mask_bs = (xx - N/2).^2 + (yy - N/2).^2 < r_bs^2; I_meas(mask_bs) = 0;

poissrnd的输入是期望光子数，输出是带泊松涨落的随机计数值。先把强度归一化到峰值photon_count，加完噪声再缩放回去，是为了让噪声后的强度仍处于原仿真量级，既不过度溢出也不丢失动态范围。中心遮挡半径 3 个像素，对应实验中阻挡直射光的小圆盘。遮挡区域后续在迭代里要特殊处理，不能把它当成真正的零强度测量值。

3.3 自相关支撑初始化：从强度图猜物体位置

相位恢复不能从空支撑出发，有个经典做法是通过强度图的自相关来估计支撑位置和尺寸。自相关的支撑宽度大约是物体支撑的两倍，所以先对测量振幅做一次逆傅里叶变换，取显著区域再缩小一半。

% 自相关，得到物体支撑的粗估计 ac = fftshift(ifft2(ifftshift(sqrt(I_meas)))); ac = abs(ac); th = 0.25 * max(ac(:)); blob = ac > th; % 显著区域，约等于物体与其翻转的卷积 % 区域质心作为物体中心，等效边长折半作为物体尺寸估计 stats = regionprops(blob, 'Centroid', 'Area'); cx = round(stats.Centroid(2)); cy = round(stats.Centroid(1)); L_ac = sqrt(stats.Area); % 自相关等效边长 obj_d = max(6, round(L_ac / 2)); % 物体半径估计 % 生成圆形初始支撑 support = false(N); [sx, sy] = meshgrid(1:N, 1:N); support((sx - cx).^2 + (sy - cy).^2 < obj_d^2) = true;

regionprops来自图像处理工具箱，如果环境中没有这个工具箱，可以用mean和std找质心，也可以固定取矩阵中心作为支撑中心——模拟里物体通常就在中心附近。自相关撑起的区域是物体支撑的自卷积，所以L_ac / 2是物体尺寸的合理估计。这个初始支撑不必很准，后面 HIO 迭代中会逐步修正。

3.4 HIO 迭代主循环与 shrinkwrap 支撑更新

核心迭代采用混合输入输出（HIO）算法，配合每隔一定轮数更新一次支撑的 shrinkwrap 策略。HIO 在支撑外的处理比误差还原（ER）多了一个负反馈项，能有效摆脱局部极小，是 CDI 模拟源码里最常见的主循环。

beta = 0.8; % HIO 反馈参数 n_iter = 500; rng(1); obj_est = sqrt(I_meas) .* exp(1i * 2 * pi * rand(N)); obj_est = fftshift(ifft2(ifftshift(obj_est))); % 随机相位初始 err_f = zeros(n_iter, 1); h = fspecial('gaussian', [5 5], 1.0); % shrinkwrap 平滑核 for k = 1:n_iter % 傅里叶域：保持测量振幅，替换相位 F_est = fftshift(fft2(obj_est)); F_upd = sqrt(I_meas) .* exp(1i * angle(F_est)); g = fftshift(ifft2(ifftshift(F_upd))); % 回到实空间 % 实空间：支撑内取 g，支撑外按 HIO 规则更新 obj_new = g .* support ... + (obj_est - beta * g) .* (~support); obj_est = obj_new; % 傅里叶域相对误差 err_f(k) = norm(abs(F_est) - sqrt(I_meas), 'fro') ... / norm(sqrt(I_meas), 'fro'); % 每 20 轮更新支撑（shrinkwrap） if mod(k, 20) == 0 blurred = imfilter(abs(obj_est), h, 'replicate'); support = blurred > 0.15 * max(blurred(:)); support = imdilate(support, strel('disk', 2)); end end

这段循环里两个约束交替生效：傅里叶域用测量振幅替换估计振幅、保留估计相位；实空间则把支撑域内的值直接替换为g，支撑域外的值按obj_est - beta * g更新。beta越大，支撑外的抑制越强，收敛快但容易振荡；beta太小则摆脱局部极小的能力弱。0.8 是多数模拟任务里不用怎么调的默认值。

支撑更新用的是高斯模糊后取阈值，再做半径为 2 的膨胀，避免支撑边界收缩得过于激进。初始支撑偏大的情况下，这个策略会在迭代中逐渐把支撑拉近到物体真实边界。

注意：中心遮挡mask_bs区域的强度是人为置零的，不是真实测量。严格处理时该区不应参与傅里叶约束，更稳的做法是给这块区域权重 0。简单模拟里直接让sqrt(I_meas)为零也能收敛，但重建物体中心会有一个轻微暗斑，这是 beamstop 伪影而不是算法问题。

4. 相位恢复算法的参数设定、收敛判据与抗噪处理

4.1 ER、HIO、RAAR 三种更新的差别

在实空间投影这一步，不同算法的差异只有一两行代码，收敛行为却完全不同。误差还原（ER）是支撑外直接置零，单调去逼近一组可行解，但很容易在第一个局部极小值附近停住。HIO 引入了反馈项，让支撑外的残余误差反向作用到下一次估计上，跳出局部极小的能力明显更强。RAAR 则在 HIO 与反射型更新之间做插值，对噪声和高饱和度区域更稳健。

MATLAB 里从 HIO 换到 ER 只是把支撑外的一行换成obj_est = obj_est .* support，但实际使用时不要单跑 ER，常见做法是先跑 100 轮 HIO 让支撑收敛，再切到 ER 做最后的平滑细化。RAAR 实现更复杂一些，等位相恢复遇到明显噪声时再考虑替换，入门阶段把 HIO 调好就够用。

4.2 beta、迭代轮数与收敛判据的配置

beta 不是越大越好，也不是越小越稳。从经验看，beta 在 0.7 附近时 HIO 的振荡幅度适中；降到 0.5 以下，每次迭代对支撑外误差的修正太小，500 轮下来误差降得又慢又不彻底；升到 1.0 以上，误差曲线会出现周期性震荡，看起来像在两组解之间反复横跳。

迭代轮数的判断要结合误差曲线。傅里叶域相对误差err_f的表达式已经在第 3 章代码里定义，正常的收敛过程是前 50 轮快速下降，然后进入平缓滑行。如果发现曲线走平后还要继续加轮数去“碰运气”，正确做法是换一个随机相位初值重启，或者先检查支撑是否过小。

% 收敛判定：最近100轮相对变化小于0.5%就提前停下 if k > 100 recent = err_f(k-99:k); if (max(recent) - min(recent)) / max(recent) < 5e-3 fprintf('在第 %d 轮达到收敛\n', k); break; end end

这个判据只适用于模拟里固定真值的场景。真实实验没有真值可对比，只能靠傅里叶域误差走平来判断迭代是否进入稳态。如果连走平都看不到，优先检查过采样比和支撑初始化，而不是调 beta。

4.3 噪声、饱和像素与权重掩膜的抗噪策略

泊松噪声是 CDI 模拟里最接近实验实际的噪声模型。一个值得记住的结论是：在sqrt(I)域做傅里叶约束，而不是在I域做，本身就接近泊松噪声的极大似然处理，这也是主循环里都写sqrt(I_meas)的原因，不是图方便，而是有统计依据。

对于饱和像素和 beamstop 区域，处理思路应该彻底反过来：完全不约束，让相位自由演化。实现方法是用一个权重掩膜控制傅里叶域“哪些像素被强迫等于测量振幅”，权重为 0 的位置跳过约束。

weight = ones(N); weight(mask_bs) = 0; % 中心遮挡区不约束 sat_pix = I_meas > 0.99 * max(I_meas(:)); % 饱和区也不约束 weight(sat_pix) = 0; for k = 1:n_iter F_est = fftshift(fft2(obj_est)); amp_target = sqrt(I_meas) .* weight ... + abs(F_est) .* (1 - weight); F_upd = amp_target .* exp(1i * angle(F_est)); % 后续实空间更新与第3章相同 end

这段代码把“振幅替换”改为“振幅混合”，权重为 0 的像素保留当前估计振幅，而不是强行拉向测量值。加了权重掩膜后，beamstop 伪影和饱和点扩散都会明显减弱，代价是有效约束像素变少，迭代需要更多轮数。

5. 重建质量怎么量化：误差指标、伪影识别与 shrinkwrap 自适应

5.1 模拟场景下可用的三个硬指标

模拟的好处是有真值，可以直接算误差。三个指标里，归一化 RMSE 反映振幅的整体偏差，相关系数反映结构相似度，傅里叶域误差则是迭代收敛性的直接依据。代码很短：

mask = support & (abs(obj) > 0.1 * max(abs(obj(:)))); rmse = sqrt(mean(abs(obj(mask) - obj_est(mask)).^2 ./ abs(obj(mask)).^2)); cc_amp = corr2(abs(obj), abs(obj_est));

corr2返回 -1 到 1 的值，0.98 以上可以算高质量重建；RMSE 需要看具体物体动态范围，一般小于 0.1 说明振幅恢复得不错。值得注意的是，只用angle(obj_est)直接对比相位没有意义，因为整体相位还存在一个任意常数偏移，对比前要先把重建相位减去平均值。

5.2 从伪影形态反推参数问题

重建结果出现两类常见伪影时，别急着换算法。第一类是背景散斑噪声，表现为支撑外有零散亮点，通常是初始支撑过大，迭代中 shrinkwrap 阈值太低没有及时收缩支撑。第二类是物体边缘振铃和条纹，通常是支撑过小，物体真实边界被截断了，重建被迫把能量挤进了一个过小的框。

一个快速判断技巧是显示重建振幅的动态范围。支撑过大时，最大振幅会明显超出真值；支撑过小时，物体内部会出现波纹状明暗交替，且边缘处尤其明显。对应修法是调整 shrinkwrap 的阈值参数，而不是重跑整个模拟。

5.3 shrinkwrap 阈值与膨胀半径的调节技巧

shrinkwrap 实际只有三个可调量：高斯平滑核尺寸、阈值、膨胀半径。核太大会把支撑边界抹得模糊，核太小则容易把目标内部噪声识别成支撑区域。推荐从[5 5]核、标准差 1.0、阈值 0.15、膨胀半径 2 起步，这四个值覆盖了多数模拟场景。噪声明显偏高时，把阈值从 0.15 提到 0.2，同时把膨胀半径从 2 增加到 4，撑起一段缓冲带，支撑会更稳。实际跑模拟时，我一般固定 beta 为 0.8，观察误差曲线在 200 轮后的走势，再决定优先调阈值还是调膨胀半径，这比盲目增加迭代轮数有效得多。

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