基于DMD实现透过多模光纤(MMF)的聚焦

前言

在波前整形(WFS)领域，虽然LC-SLM能实现聚焦PBR系数 $\frac{\pi }{4}$ 的相位调制，但其典型值仅为 $100$ Hz的低刷新率，导致无论是基于迭代的波前优化还是测TM实现透过复杂介质聚焦均较慢，从而限制了其应用场景。 相比之下，含大量数字微镜阵列的DMD具备高达23KHz的刷新率，其每个微镜都能在“ON”和“OFF”两种状态间切换，能提供聚焦PBR系数 $\frac{1}{2\pi }$ 的二元幅度调制。基于DMD的波前优化算法[1,2] 和传输矩阵™测量用于透过光学散射实现聚焦[3,4,5] 已经有不少研究。这篇博客总结一下基于DMD实现透过多模光纤(MMF)的聚焦。

基于DMD表征MMF的背景研究

同散射介质类似，MMF也是一种复杂介质。光透过MMF传播会发生模式色散因而产生波前畸变，导致其无法直接用于聚焦，或形成object图案用于成像。通常基于全息法让参考光与MMF输出散斑干涉来还原光场相位信息，进而测量MMF的复值TM，但全息法装置复杂，且对系统稳定性依赖高。对于基于DMD表征MMF的非全息法，有：

1. 复值TM：基于变分法的相位复原(phase retrieval)，仅根据MMF输出光场的强度测量和DMD二值化输入图案估计复值TM[3]。
2. 部分TM：包括①对基于Hadamard基的DMD pattern巧妙构造，仅由输出光场强度直接测出复值TM的实部，得到实值TM[6,10]；②基于一系列输入掩膜下输出光场强度的变化测量二值TM，用于快速的DMD二元幅度调制[7]。
3. DMD微镜开闭的条件概率算法(CPA)[8]，用于决定最优DMD pattern实现聚焦。
4. 基于DMD的迭代波前优化算法，典型的有遗传算法 (GA)[2]、粒子群优化(PSO)算法、蚁群(AC)算法以及针对波前介质扰动改进后的动态优化算法(DMA)[9]等

针对上述四类测量TM或者迭代优化波前的方法，文献[11]比较了各自的理论聚焦PBR，以及在相同实验装置和条件下产生最优DMD pattern进行聚焦的速度（仿真和实验），同时本文也主要是基于翻译该文献加上个人总结而得来。

测复值TM的理论聚焦PBR

一旦测出复值TM，要在MMF输出光场第 $m$ 个mode处聚焦，我们就取TM第 $m$ 行元素 $t_{mn}(n=1,2,...N)$ 二值化处理构成DMD pattern，有两种方式：

1. 相位在 $(-\frac{\pi }{2},\frac{\pi }{2})$ 之间（实部>0）的对应1，否则对应0；
2. 相位在 $(0,\pi )$ 之间（虚部>0）的对应1，否则对应0；

这样DMD输入光场能透过MMF在目标输出位置发生干涉相长形成聚焦。事实上，打开某些微镜对背景光强的比对目标光强的贡献还大，因此，上述两个决定最优DMD二值图案的标准(criteria) 需要修正，当关闭这些不该打开的微镜能进一步提升聚焦PBR，图解如下：

令 $N$ 为DMD微镜 (input modes) 数，假设 $t_{mn}$ 服从Rayleigh分布且 ${\left|t_{mn}\right|}^2$ 遵循指数衰减 ${\left|t_{mn}\right|}^2\sim e^{-1/2{\sigma }^2}$ ，其中 $2{\sigma }^2$ 是$t_{mn}$每个元素的总体平均强度，$\phi$ 为绝对相位差异的上确界：$0⩽\left|\Delta \theta \right|<\phi$，则基于DMD二元幅度调制的理论聚焦PBR[12]为： $$PBR\approx \frac{(N-1)si{n}^2\phi }{4\phi }$$
其中绝对相位差异的上确界$\phi =\pi /2$。注意到该理论PBR公式只考虑 $t_{mn}$ 的相位信息用于决定最优DMD pattern做聚焦，适用于测复值TM；其它测实值TM、条件概率算法CPA以及GA等迭代法能实现的理论聚焦PBR将会与之做对比。

测实值TM的改进算法

测实值TM是由Yongkeun Park在2017年最先提出的[6]，这里先简要归纳如下：

利用Hadamard基 ${\mathbf{h}}_{\mathbf{n}}$ 作为DMD输入，对应的输出光场为 ${\sigma }_{\mathbf{n}}$，输入与输出的线性关系可由复值TM $\mathbf{T}$表示为：
$$\left[{\sigma }_1,{\sigma }_2\cdots {\sigma }_{\mathbf{N}}\right]=\mathbf{T}\left[{\mathbf{h}}_1,{\mathbf{h}}_2\cdots {\mathbf{h}}_{\mathbf{N}}\right]$$
复值TM实部 $Re(\mathbf{T})=\mathbf{S}{\mathbf{H}}^T=\mathbf{S}\mathbf{H}$，其中 $\mathbf{S}=\frac{1}{N}Re\left[{\sigma }_1,{\sigma }_2\cdots {\sigma }_{\mathbf{n}}\right]$，那么怎么获得实值TM $Re(\mathbf{T})$呢？
对于Hadamard矩阵 $\mathbf{H}=\left[{\mathbf{h}}_1,{\mathbf{h}}_2\cdots {\mathbf{h}}_{\mathbf{n}}\right]$，我们构造如下两个二值化Hadamard基作为DMD pattern：
$${\mathbf{v}}_{\mathbf{n}}^{\pm }=\frac{1}{2}({\mathbf{h}}_1\pm {\mathbf{h}}_{\mathbf{n}}),n=1,2\cdots N$$
对应 ${\mathbf{v}}_{\mathbf{n}}^{\pm }$ 的输出光强 ${\mathbf{I}}_{\mathbf{n}}^{\pm }$ 为：
$${\mathbf{I}}_{\mathbf{n}}^{\pm }={\left|\mathbf{T}{\mathbf{v}}_{\mathbf{n}}^{\pm }\right|}^2={\left|\frac{1}{2}\mathbf{T}({\mathbf{h}}_1\pm {\mathbf{h}}_{\mathbf{n}})\right|}^2=\frac{1}{4}{\left|{\sigma }_1\right|}^2+\frac{1}{4}{\left|{\sigma }_{\mathbf{n}}\right|}^2\pm \frac{1}{2}Re({\sigma }_1^{\ast }{\sigma }_{\mathbf{n}})(复数相加的模的平方运算)$$
令 ${\sigma }_1$ 相位为 0 不失一般性，则我们得到输出光场 ${\sigma }_{\mathbf{n}}$ 的实部为 $Re\left[{\sigma }_{\mathbf{n}}\right]=\frac{{\mathbf{I}}_{\mathbf{n}}^{+}-{\mathbf{I}}_{\mathbf{n}}^{-}}{\sqrt{{\mathbf{I}}_1}}$ ，其中 ${\mathbf{I}}_1={\left|\mathbf{T}{\mathbf{v}}_1\right|}^2={\sigma }_1^2$
这样我们就能通过组合得到二值化的Hadamard基，由输出光场强度直接测出复杂介质的实值TM，要实现透过复杂介质在第 $m$ 个mode处聚焦，取实值TM第 $m$ 行，按照元素>0的取1否则取0的方法得到最优DMD输入图案。

Tianrui ZHAO等人则进一步提出了实值强度TM ($RVITM$)，连接二值DMD输入与输出光强[10]。

同Park测实数值TM构造二值化Hadamard基一样，构造 两个大小为$N\times N$二值化Hadamard矩阵${\mathbf{H}}_1=(\mathbf{H}+1)/2,{\mathbf{H}}_2=(-\mathbf{H}+1)/2$，其中每个列向量作为Hadamard基构成DMD pattern，则二值输入与输出光强的线性关系为：
$$\left[\begin{array}{cccc}& I_1^1& \cdots & I_1^{2N}\\ & ⋮& \ddots & ⋮\\ & I_M^1& \cdots & I_M^{2N}\end{array}\right]=RVITM\cdot \left[{\mathbf{H}}_1,{\mathbf{H}}_2\right]$$
先对上式都右乘主对角元为2的对角矩阵，再考虑到 ${\mathbf{I}}_1=RVITM\cdot {\mathbf{h}}_1$， 其中 ${\mathbf{h}}_1$ 为全1向量，若等式右边矩阵$\left[2{\mathbf{H}}_1,2{\mathbf{H}}_2\right]$ 减去1，则有：
$$\left[\begin{array}{cccc}& 2I_1^1-I_1^1& \cdots & 2I_1^{2N}-I_1^1\\ & ⋮& \ddots & ⋮\\ & 2I_M^1-I_M^1& \cdots & 2I_M^{2N}-I_M^1\end{array}\right]=RVITM\cdot \left[2{\mathbf{H}}_1-1,2{\mathbf{H}}_2-1\right]=RVITM\cdot \left[\mathbf{H},-\mathbf{H}\right]$$
可见等式右侧变成 $RVITM$ 直接与大小为$2N\times N$ 的Hadamard矩阵 $\left[\mathbf{H},-\mathbf{H}\right]$相乘！因此我们能对二值化Hadamard矩阵 ${\mathbf{H}}_1$ 和 ${\mathbf{H}}_2$ 作为输入的MMF输出光强做一些处理，由Hadamard矩阵性质：${\left[\mathbf{H},-\mathbf{H}\right]}^{-1}=\frac{1}{2N}\cdot {\left[\mathbf{H},-\mathbf{H}\right]}^T$，我们可以直接算出实值强度TM为：
$$RVITM=\frac{1}{2N}\cdot \left[\begin{array}{cccc}& 2I_1^1-I_1^1& \cdots & 2I_1^{2N}-I_1^1\\ & ⋮& \ddots & ⋮\\ & 2I_M^1-I_M^1& \cdots & 2I_M^{2N}-I_M^1\end{array}\right]\cdot {\left[\mathbf{H},-\mathbf{H}\right]}^T$$
进一步，通过对测出的$RVITM$求伪逆，能用于从MMF输出散斑中重建输入图像：
$$I_{image}=RVIT{M}^{+}\cdot I_{out}\approx I_{in}$$

可以发现：实值强度TM ($RVITM$)较实值TM的优势是直接连接二值DMD输入与输出光场的强度测量值，而非输出光场的实部，从而能直接从MMF输出散斑测量结果中重建输入图像。

Fig 2. (a) 用RVITM表征MMF图解, (b)基于RVITM求逆从MMF输出散斑图案重建输入图像(KCL), （c) 测RVITM的实验设置

进一步， 实值强度TM ($RVITM$) 还能表达成编码相应复值TM的幅度与相位信息的形式[11].

由于输出光场强度与其激发的PA信号信号幅度有 系数为$\alpha$ 的正比例关系，因此考虑测量结果为目标处输出光强时，连接第 $n$ 个DMD微镜处输入光场与 $m$ 处目标输出光强的 $RVITM$元素为：
$$rvit{m}_{mn}=A_{mn}{A}_{R}cos({\theta }_{mn}-{\varphi }_R)$$

根据上述推导，当 ${\theta }_{mn}$ 在$\left[{\varphi }_R-\pi /2,{\varphi }_R+\pi /2\right]$ 之间，则 $rvit{m}_{mn}⩾0$，打开相对应的DMD微镜则在目标输出处发生干涉叠加用于聚焦。作者假设对应更大 $rvit{m}_{mn}$值的DMD微镜对焦点处光强比背景光强的贡献更大，因此要优先打开。作者因而按照对应 $rvit{m}_{mn}$ 大小对DMD微镜做了降序排列，在仿真和实验上探究只打开前 $P\%$ 的微镜对聚焦PBR的影响，结果如下：

基于变分法的相位复原估计复值TM

条件概率算法(CPA)

迭代波前优化 (以GA为例)

参考文献

[1] D. Akbulut, T. J. Huisman, E. G. van Putten, W. L. Vos, and A. P. Mosk, “Focusing light through random photonic media by binary amplitude modulation,” Opt. express 19, 4017–4029 (2011).
[2] X. Zhang and P. Kner, “Binary wavefront optimization using a genetic algorithm,” J. Opt. 16, 125704 (2014).
[3] A. Drémeau, A. Liutkus, D. Martina, O. Katz, C. Schülke, F. Krzakala, S. Gigan, and L. Daudet, “Reference-less measurement of the transmission matrix of a highly scattering material using a dmd and phase retrieval techniques,” Opt. express 23, 11898–11911 (2015).
[4] G. Huang, D. Wu, J. Luo, Y. Huang, and Y. Shen, “Retrieving the optical transmission matrix of a multimode fiber using the extended kalman filter,” Opt. Express 28, 9487–9500 (2020).
[5] M. N’Gom, T. B. Norris, E. Michielssen, and R. R. Nadakuditi, “Mode control in a multimode fiber through acquiring its transmission matrix from a reference-less optical system,” Opt. letters 43, 419–422 (2018).
[6] H. Yu, K. Lee, and Y. Park, “Ultrahigh enhancement of light focusing through disordered media controlled by mega-pixel modes,” Opt. express 25, 8036–8047 (2017).
[7] X. Tao, D. Bodington, M. Reinig, and J. Kubby, “High-speed scanning interferometric focusing by fast measurement of binary transmission matrix for channel demixing,” Opt. express 23, 14168–14187 (2015).
[8] T. Zhao, L. Deng, W. Wang, D. S. Elson, and L. Su, “Bayes’ theorem-based binary algorithm for fast reference-less calibration of a multimode fiber,” Opt. express 26, 20368–20378 (2018).
[9] Huanhao Li, etc, “Adaptive optical focusing through perturbed scattering media with a dynamic mutation algorithm,” Photon. Res. 9, 202-212 (2021)
[10] T. Zhao, S. Ourselin, T. Vercauteren, and W. Xia, “Seeing through multimode fibers with real-valued intensity transmission matrices,” Opt. Express 28, 20978–20991 (2020).
[11]. TIANRUI ZHAO, etc. Focusing light through multimode fibres using a digital micromirror device: a comparison study of non-holographic approaches(in press)
[12] D. Wang, E. H. Zhou, J. Brake, H. Ruan, M. Jang, and C. Yang, “Focusing through dynamic tissue with millisecond digital optical phase conjugation,” Optica 2, 728–735 (2015).