\providecommand{\main}{..}

\documentclass[main.tex]{subfiles}
\graphicspath{ {../images/} }
\begin{document}
% Sensorless AO section
\section{Refocusing using sensorless adaptive optics}

% Subs 1
\subsection{Introduction}

Adaptive optics (AO) is a technology used in optics to correct for aberrations which might arise in the optical system. It has been initially developed for astronomy applications, in order to correct for aberrations introduced by atmospheric turbulence \cite{sao_ao_astronomy}. While telescopes, such as Hubble and James Webb, can be sent into space, where the problem of imaging through the atmosphere is avoided, there are limitations, mainly in size and weight. Having larger telescopes is beneficial, as the resolution of imaging is directly related to the size of the capturing mirror/antenna. Atmospheric turbulence can be minimised by carefully choosing the location of such telescopes, for example on top of a mountain such as the Mauna Kea Observatories in Hawaii. However, aberrations can not be entirely avoided, some even arising from wind generating vibrations in the telescope structure.

Over time, adaptive optics has also found its use in high resolution microscopic imaging \cite{sao_ao_highres_imaging}. Optical aberrations in microscopy mainly arise due to specimen structure, due to changes in refractive index distribution, which hinder the formation of a high quality focus, which in the best case is limited by diffraction, a target of high resolution imaging. This non-uniform refractive index leads to changes in the optical wavefront as it passes through the sample.  The fundamental principle upon which AO works is modulation of the wavefront in order to compensate for any aberrations which might arise. If the inverse effect of the aberrations is applied onto the optical field, then it cancels out with the sample aberrations leading to a high-quality focus. In practice, such wavefront modulations are achieved using devices such as liquid crystal spatial light modulators (SLM), deformable mirrors or digital micromirror devices (DMD).

The second main component of AO is a method for characterising the aberrations, such that the correct modulation is applied.  Sensor based AO relies on the use of a wave-front sensor, which measures the amount and spatial distribution of aberrations, information which is then used to infer the required modulation for correction. One example of such sensor is the Shack-Hartmann wave-front sensor \cite{sao_shack_hartman}, which uses an array of lenslets to sample the optical field by focusing small sections onto a camera device. The displacement of the focus from the centre is then proportional to the phase gradient of the optical field, which can then be approximated using a set of local tip/tilts. Tip and tilt are linear phase gradients applied to the optical field in perpendicular transversal directions. The name arises from the possibility of inducing such changes by tipping or tilting a mirror on the optical path.

While such methods are able to directly describe the nature of the aberrations, they typically involve an increase in complexity of the optical system. The more stringent requirement is access to the aberrated optical field, which needs to be observed by the sensor. In the \textit{in vivo} imaging case, it might not be possible, as the aberrated fields occur inside the specimen.

An alternative is sensorless based AO, which uses the image acquired to infer the state of the optical field. The method functions by measuring a chosen parameter, such as intensity or contrast, and modulating the input light field, the chosen parameter is optimised.

The modulation to the field is usually represented in a basis which suits the aberrations intended to be corrected. One such basis is Zernike polynomials \cite{sao_zernike}, which describe usual aberrations such as tilt, defocus, astigmatism, etc. However, in the case of the bent fibre, the aberrations introduced do not come in the usual forms described by Zernike polynomials. As the modulation applied is phase-only in nature and sampled in a pixelated format, a basis which is defined using these constraints is more appropriate, such as the basis defined by Walsh-Hadamard modes \cite{sao_lattice_geometry}.

In this chapter, recovering the focus at a reflector at the distal facet of the optical fibre was investigated. This was done in preparation for using the memory effect to achieve correction over the entire distal facet and later over the entire target focal plane, as described in the following two chapters.

% Subs 2
\subsection{Mathematical description of Walsh mode modulation}

Correcting for the change in light transmission through a bent MMF requires the observation and description of the optical field. Although in its raw representation, the field is sampled on pixels measured by a camera, it is usually best to represent the field using a modal expansion, as explained above. However, such Cartesian grid is inherently two-dimensional. In this configuration, any linear transformation on the optical field, such as phase modulation, would be 4-dimensional in nature, as a mapping from one 2D space to another 2D space. The easiest simplification of the problem is to map the 2D arrangement of pixels onto a 1D vector. For example, this can be achieved by restructuring the data in a column only, by concatenating consecutive columns from the grid vertically. The advantage of using the mapping is that linear transformations can now be represented by a 2D matrix, which is easier to work with. This representation is not unique, any permutation of the pixel values is usable. However, whichever representation is used, it needs to be adhered to for any calculations going forward.

The field itself is of complex value. However, the modulation device (SLM) only changes the phase parameter. This is modelled by multiplying the original field elements with phase-only complex values, as illustrated in equation \ref{eq:phase_mod}, where $x_c$ is the corrected field, $x_1$,$x_2$,…,$x_p$ are the original field elements and $e^{i\theta_1}$,$e^{i\theta_2}$,…,$e^{i\theta_p}$ are the modulations applied to each pixel. The element-wise multiplication can be represented as a linear transformation of the field, represented as a column, by matrix multiplication with a diagonal matrix representation of the modulation. Conversely, as the object of interest is the modulation and not the previously calibrated field, which is a starting constant during the correction procedure, the diagonal matrix and column representation can be switched in the form below.

\begin{equation} \label{eq:phase_mod}
    x_c = 
    \begin{pmatrix}
        x_1 \cdot e^{i\theta_1} \\
        x_2 \cdot e^{i\theta_2} \\
        \vdots \\
        x_p \cdot e^{i\theta_p}
    \end{pmatrix}
    =
    \begin{pmatrix}
        x_1 & 0 & \cdots & 0 \\
        0 & x_2 & \cdots & 0 \\
        \vdots & \vdots & \ddots & \vdots \\
        0 & 0 & \cdots & x_p
    \end{pmatrix}
    \begin{pmatrix}
        e^{i\theta_1} \\
        e^{i\theta_2} \\
        \vdots \\
        e^{i\theta_p}
    \end{pmatrix}
\end{equation}

The optical field is in itself a highly dimensional set of data, each pixel representing a dimension. One method of characterising the change of such field on the output of the optical system, is to modulate one dimension at a time. However, modulating a single pixel has minimal effect on the intensity at the output. Thus, it is better to use a modal representation in which all input pixels are modulated in every dimension. One such basis transformation can be done with the use of Walsh-Hadamard transform \cite{sao_lattice_geometry}. The optical field modulation can be represented as a superposition of Walsh modes, which in terms have a phase coefficient associated.

The Walsh modes are a set of orthogonal basis functions which take only binary values of +1 and -1, on sub-intervals defined by dyadic fractions (the denominator is a power of 2). In other words, in the 2-dimensional representation, if a row/column of one Walsh function is considered, the length of each consecutive sequence of identical valued pixels must be a dyadic fraction of the length of the row/column. For example, in the $8\times8$ case shown in figure 4, all bands have a vertical/horizontal length of 1, 2, 4 or 8 pixels, associated with 1/8, 1/4, 1/2 and 1/1 dyadic fractions respectively.
One of the limitations that Walsh mode representation introduces is the size of the optical field. As the modes are generated using dyadic functions, the size of the field in pixels, must be a power of 2, such as $8\times8$ shown in figure \ref{fig:all_8x8_walsh}.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{images/walsh_exmaples_8x8.png}
    \caption{All the $8\times8$ size Walsh modes. Black represents a value of -1 and brown represents a value of +1.}
    \label{fig:all_8x8_walsh}
\end{figure}

Through the mapping between a 2D matrix and a column vector, there is an equivalence between the inherent 1D representation of Walsh modes and the 2D representation in Fig. \ref{fig:all_8x8_walsh}. The 1D column representation can be stacked into a square matrix, leading to the generation of the Walsh-Hadamard matrix, usually denoted with $\mathbf{H}$. Multiplying it with a column vector containing the coefficients, the phase field is generated in the pixel domain is calculated. This is represented in equation \ref{eq:hadamard}, where $\theta_1$,…,$\theta_p$ are the phase values of the optical field pixels, $\phi_1$,…,$\phi_W$ are the Walsh mode coefficients attributed to each Walsh mode, $p$ is the number of pixels and $W$ is the number of Walsh modes $(p=W)$.

\begin{equation} \label{eq:hadamard}
    \begin{pmatrix}
        \theta_1 \\
        \theta_2 \\
        \vdots \\
        \theta_p
    \end{pmatrix}
    =
    \mathbf{H}
    \begin{pmatrix}
        \phi_1 \\
        \phi_2 \\
        \vdots \\
        \phi_W
    \end{pmatrix}
\end{equation}

It is important to note that the Walsh transformation is applied to the numerical phase values of the optical field ($\theta_1$,$\theta_2$,…), while the modulation to the field is done using the complex equivalent ($e^{i\theta_1}$, $e^{i\theta_2}$,…).

The main advantage of using Walsh modes basis instead of the pixel basis, is that changing one coefficient is equivalent to a phase change across the entire optical field. Thus, changes in the output intensities are larger and measurement is easier to perform. It is important to note that the Walsh-Hadamard basis is by no means the only basis which has the properties of modulating a large area of the input field, however, due to the basis generation being readily available in MatLab, it has been chosen to be investigated. 

In the next steps, the Walsh mode representation is applied to the optical field transmitted through a MMF. It will be shown that changing only one Walsh mode coefficient leads to a sinusoidal evolution of the intensity of a chosen output pixel. This is not a property specific to the Walsh modes, but rather to any modulation between two arbitrary wavefronts, where only the global phase of one of the wavefronts is changed. This property will be useful in creating an optimisation algorithm to correct the change of transmission through a bent fibre. Introducing the modulation, as described in equation \ref{eq:phase_mod}, into the TM representation leads to the following equation, where $\mathbf{y}$ is the output field, $\mathbf{T}$ is the transmission matrix between the modulated field and the output field and $\mathbf{x_0}$ is the unmodulated input field:

\begin{equation} \label{eq:mod_to_TM}
    \mathbf{y} = \mathbf{T} \mathrm{diag}(\mathbf{x_0}) \exp{i \mathbf{\theta}} = \mathbf{T_{x0}} \exp{i \mathbf{\theta}}
\end{equation}

This is the general form of a phase-only modulation for the transmission of light through a MMF. This can alternatively be expressed using a Walsh mode representation, as defined using the form in equation \ref{eq:hadamard}, where $\phi_1$,$\phi_2$,… are the phase coefficients of the Walsh modes and \textbf{H} is the Walsh-Hadamard matrix.

As the target of the optimisation is a specific output, the next step is to focus on one output pixel (of index $s$), and focus equation \ref{eq:phase_mod} to the specific row of \textbf{T}, associated with the chosen output. In the equations below, $y_s$ is the output at pixel $s$, $(t_{x0})_{sk}$ is the element at position $(s,k)$ in the matrix $\mathbf{T} \mathrm{diag}(\mathbf{x_0})$ and $(h)_{kq}$ is the element at position $(k,q)$ in \textbf{H}, the Walsh-Hadamard matrix. $P$ is the number of pixels and $W$ is the number of Walsh modes. While $P$ and $W$ are equal, they are still noted with different letters in order to differentiate between summation over pixels and summation over Walsh modes.

\begin{equation} \label{eq:Walsh_element_first}
    y_s = \sum_{k=1}^P(t_{x0})_{sk} \cdot e^{i\theta_k}
\end{equation}

\begin{equation} \label{eq:Walsh_element_second}
    \theta_k = \sum_{q=1}^W (h)_{hq} \phi_q
\end{equation}

\begin{equation} \label{eq:Walsh_element_third}
    y_s = \sum_{k=1}^P (t_{x0})_{sk} \cdot \exp \left[ i \left( \sum_{q=1}^W (h)_{kq} \cdot \phi_q \right) \right]
\end{equation}

Combining the TM expression for a specific output in eq. \ref{eq:Walsh_element_first} with the sum expression of the matrix multiplication in eq. \ref{eq:hadamard}, reduced to a specific pixel as given in eq. \ref{eq:Walsh_element_second}, results in the form in eq. \ref{eq:Walsh_element_third} which expresses a specific output in terms of the Walsh coefficients which form the modulation.

As described previously, the point of interest of Walsh mode optimisation is the evolution of the output pixel as a function of only one Walsh mode coefficient. Hence, the desired Walsh mode coefficient (of index $v$, coefficient $\phi_v$) is extracted out of the summation in eq. \ref{eq:Walsh_element_third} as the other terms are independent of this coefficient.

\begin{equation} \label{eq:Walsh_element_fourth}
    y_s = \sum_{k=1}^P (t_{x0})_{sk} \cdot \exp \left[ i \left( (h)_{kv} \cdot \phi_v + \sum_{q=1, q\neq v}^W (h)_{kq} \cdot \phi_q \right) \right]
\end{equation}

\begin{equation} \label{eq:Walsh_element_fifth}
    y_s = \sum_{k=1}^P (t_{x0})_{sk} \cdot \exp \left[ i \left( \sum_{q=1, q\neq v}^W (h)_{kq} \cdot \phi_q \right) \right] \cdot \exp [ i \cdot (h)_{kv} \cdot \phi_v ]
\end{equation}

Addition in the exponential function in eq. \ref{eq:Walsh_element_fourth} is expanded to multiplication of exponential functions in eq. \ref{eq:Walsh_element_fifth}. The first two elements of the multiplication in eq. \ref{eq:Walsh_element_fifth} are constant with respect to Walsh coefficient $\phi_v$. Thus, for simplicity, they can be collapsed into a general term $z_k$.

\begin{equation}
    y_s = \sum_{k=1}^P z_k \cdot \exp [ i \cdot (h)_{kv} \cdot \phi_v ]
\end{equation}

The next step is to make use of the property of a Walsh mode, that it contains a set of $+1$ and a set of $-1$ $(h)$ terms. 
These can be separated into two summations (set of $k_+$ and set of $k_-$ indices, which indicate the +1 and -1 terms, respectively), as it allows for the $(h)_{kv}$ terms to be substituted with $+1$ or $-1$, allowing for the exponent terms to be taken out of the summation, the final summations over $k$ containing only the constant terms described before. The result in eq. \ref{eq:Walsh_element_sixth} is only the sum of two terms, where the only variable is $\phi_v$:

\begin{gather}
    y_s = \sum_{k_+} z_{k_+} \cdot \exp [i \cdot (+1) \cdot \phi_v] + \sum_{k_-} z_{k_-} \cdot \exp [i \cdot (-1) \cdot \phi_v] \nonumber \\
    y_s = \exp (i \cdot \phi_v) \cdot \sum_{k_+} z_{k_+} + \exp (-i \cdot \phi_v) \cdot \sum_{k_-} z_{k_-} \nonumber \\
    y_s = \exp (i \cdot \phi_v) \cdot z_+ + \exp (-i \cdot \phi_v) \cdot z_- \label{eq:Walsh_element_sixth}
\end{gather}

The next steps to follow have the aim of proving the sinusoidal evolution of the output pixel intensity as a function of Walsh mode coefficient $\phi_v$. The expression in equation \ref{eq:Walsh_element_sixth} describes the evolution of the output pixel as a function of the coefficient of Walsh mode $v$. $z_+$,$z_-$ are two complex numbers which arise from the summation and are dependent on the TM, non-corrected input and the coefficients of the other Walsh modes. The value of interest is the intensity of output $y_s$, which is independent of the global phase. Thus, the global phase reference value can be chosen such that the phases of $z_+$  and $z_-$ are opposites in order to make the calculations simpler. 

\begin{gather}
    z_+ = |z_+| \cdot e^{i\rho} \nonumber \\
    z_- = |z_-| \cdot e^{-i\rho} \nonumber
\end{gather}

Substituting the new forms of $z_+$  and $z_-$ into eq. \ref{eq:Walsh_element_sixth} and solving for the output intensity leads to:

\begin{gather}
    y_s = \exp (i \cdot \phi_v) \cdot |z_+| \cdot e^{i\rho} + \exp (-i \cdot \phi_v) \cdot |z_-| \cdot e^{-i\rho} \nonumber \\
    y_s = (|z_+| + |z_-|) \cdot \cos (\phi_v + \rho) + i \cdot (|z_+| - |z_-|) \cdot \sin (\phi_v + \rho) \nonumber \\
    |y_s|^2 = ( |z_+|^2 + |z_-|^2 ) + 2 \cdot |z_+| \cdot |z_-| \cdot \cos [2(\phi_v + \rho)] \label{eq:Walsh_element_seventh}
\end{gather}

The equation above describes a sinusoidal evolution with fixed frequency ( 2 radians/Walsh mode coefficient ), irrespective of fibre and input parameters, but variable offset ($OF = |z_+|^2 + |z_-|^2$), amplitude ($A = 2 \cdot |z_+| \cdot |z_-|$) and phase offset ($PH = 2 \rho$).

\begin{equation}
    |y_s|^2 = OF + A \cdot \cos(2 \phi_v + PH)
\end{equation}

The sinusoidal evolution of the target intensity was verified using the simulation program and it is shown in figure \ref{fig:walsh_sine_sim} for two example chosen Walsh modes.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{images/walsh_sinusoidal_example.png}
    \caption{Simulated output intensity variation with respect to a Walsh mode coefficient. Modulation done at the proximal facet plane. The results illustrate the predicted sinusoidal variation.}
    \label{fig:walsh_sine_sim}
\end{figure}

The derivation in this section of the chapter proved the sinusoidal evolution of the intensity of a chosen output location with respect to the variation of a single Walsh mode coefficient. The modulation was done at the proximal facet plane. However, the same applies for the modulation at the SLM plane as well.

It is assumed that the SLM is in the Fourier plane of the proximal facet. Thus, the propagation between the SLM and proximal facet planes can be described through the use of a Fourier transform. In finitely sampled square 2D planes (same number of samples in both dimensions), the Fourier transform can be described through the use of a matrix transformation. Specifically, $\mathbf{FT} = \mathbf{F}_{shift} \cdot (\mathbf{DFT}_P \otimes \mathbf{DFT}_P)$, where $\mathbf{FT}$ is the matrix describing the Fourier transform, $\mathbf{F}_{shift}$ represents the matrix describing the Fourier shift function applied to a 2D sampled plane, $\mathbf{DFT}_P$ represents the discrete Fourier transform matrix of size $P \times P$ ($P$ being the number of samples along either dimension) and $\otimes$ being the Kronecker product.

The only change to the derivation in this subsection is that the matrix $\mathbf{T}_{x0}$ now is changed from $\mathbf{T}_{x0} = \mathbf{T} \cdot \mathrm{diag}(\mathbf{x}_0)$ to $\mathbf{T}_{x0} = \mathbf{T} \cdot \mathbf{FT} \cdot \mathrm{diag}(\mathbf{x}_0)$. Thus, the rest of the derivation is the same and the result is the same as well, though with different parameters ($OF$, $A$ and $PH$).

% Subs 3
\subsection{Calculating and applying the correction}

The target of the sensorless AO optimisation procedure in this case is to maximise the output intensity of a target pixel.  Having described the form of the output intensity as a function of Walsh mode coefficient in the previous section, the next step is to calculate the required Walsh mode coefficient in order to maximise the output intensity. In order to maximise the output intensity, the coefficient $\phi_v$ has to be set to $\frac{-PH}{2}$, a value which needs to be determined.

As there are three unknown terms ($OF$,$A$,$PH$), at least three measurements are required to determine them. This approach is named phase-shifting interferometry. While three measurements are the theoretical minimum number, in practice, the optimisation procedure is susceptible to noise. Thus, optimisation procedures with a larger number of measurements have been created (see chapter 6.7 in \cite{phase_shifting_interferometry}). However, introducing a larger number of measurements increases the time required for an iteration of the optimisation procedure, which leads to an overall increase in the total amount of time required for the procedure. Thus, the desired number of measurements will need to be chosen in practice based on a balance between the time required and the amount of noise. In this chapter, I limited the number of measurements to three, as the procedure is mostly intended to be used in a simulation environment where noise is not added.

In order to simplify the calculations, by working with simpler forms of the function, the three values to be measured are chosen as: $\phi_v = 0,-\frac{\pi}{4},\frac{\pi}{2}$, which will result in three measurements $I_1$,$I_2$ and $I_3$ respectively.

\begin{gather}
    I_1 = f(0) = OF + A \cdot \cos (PH) \\
    I_2 = f(-\frac{\pi}{4}) = OF + A \cdot \cos (PH - \frac{\pi}{2}) = OF + A \cdot \sin(PH) \\
    I_3 = f(\frac{\pi}{2}) = OF + A \cdot \cos (PH + \pi) = OF - A \cdot \cos(PH) 
\end{gather}

Using the three measurements the unknown variables can be calculated.

\begin{gather}
    OF = \frac{I_1 + I_3}{2} \\
    A = \sqrt{(I_1 - OF)^2 + (I_2 - OF)^2}   \\
    PH = \pm \cos^{-1} \left( \frac{I_1 - OF}{A} \right)
\end{gather}

PH has a $\pm$ in its form because ‘cos’ is an even function. The correct sign is identical to the sign of  $\sin(PH) = \frac{I_2 - OF}{A}$.

As the non-corrected case is when $\phi_v=0$, on the sinusoidal function, it is assured that the correction will result in an intensity which is higher or, in the worst case, equal to the starting intensity. The procedure is then repeated in an iterative process for a chosen order of Walsh modes.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{images/walsh_optimisation_diagram.png}
    \caption{Diagram description of the steps undertaken in one Walsh mode correction iteration.}
    \label{fig:walsh_opt_diagram}
\end{figure}

The advantage of this Walsh mode correction approach is that access to the distal facet is not required. If a reflector is placed at the target output location, maximising for the total reflected intensity is identical to maximising for the focus intensity at the reflector. The disadvantage is that the Walsh mode algorithm only corrects for one focus. However, the correction can be extended from one location, through the use of memory effect, discussed in more detail in chapter \ref{sec:mem_effect}.

%However, the order in which the Walsh modes are iterated through can affect the final result of the correction. As an exhaustive search of all permutations is unfeasible, a few permutations were chosen to illustrate such differences and give a broad look into the performance of Walsh mode correction. 

% subs
\subsection{Simulating the Walsh mode optimisation}

%subsubs
\subsubsection{Methodology}

While the Walsh modes represent an orthogonal basis for representing the phase modulation applied by the SLM, the order the Walsh modes are iterated through has an effect on the optimisation result.

The SLM in the optical system is a liquid-crystal phase-only type and has a resolution of $512\times512$ pixels. However, for the simulation smaller sizes have been used. Initially, a set of 4096 Walsh modes were used ($64\times64$), then this was later extended to a set of 16384 Walsh modes ($128\times128$). Smaller sizes of Walsh modes can be used on the SLM, as tiling the SLM can easily be done (i.e. $64\times64$ Walsh modes can be used on a $512\times512$ SLM, by expanding each pixel of the Walsh mode to a $8\times8$ size tile on the SLM).

The base order of the Walsh modes used is the natural (or Kronecker) order. This was generated with the use of the ``hadamard'' function in MatLab. While it is not specified in the description of the function, this order is usually generated through a recursive process. An illustration of the Hadamard matrix in the natural order is given in figure \ref{sfig:walsh_had_matrix}.

\begin{figure}
    \centering
    \begin{subfigure}[t]{0.4\textwidth}
        \includegraphics[width = \textwidth]{images/4x4_size_16_hadamard.png}
        \caption{Image of a $16\times16$ Hadamard matrix. Each column/row is a Walsh mode of size $4\times4$. Black represents -1 and brown represents +1.}
        \label{sfig:walsh_had_matrix}
    \end{subfigure}
    \begin{subfigure}[t]{0.4\textwidth}
        \includegraphics[width = \textwidth]{images/walsh_diag_perm_diagram.png}
        \caption{Diagram showing the diagonal permutation used in the Walsh mode optimisation simulations. The numbers represent the Walsh mode indices.}
        \label{sfig:walsh_diag_perm}
    \end{subfigure}
    \caption{The Walsh modes were generated in natural order. Three permutations were used in the simulations. The diagonal permutation is illustrated here.}
    \label{fig:walsh_diagrams}
\end{figure}

Three permutations were simulated: the original order of the modes, a ``diagonal'' permutation and a random permutation. 
The diagonal permutation was generated by iterating along the main diagonals of the matrix of indices, as shown in figure \ref{sfig:walsh_diag_perm}. The motive for choosing the permutation is in order to try to limit any decrease in optimisation efficiency due to similarities which might arise between Walsh modes which are either consecutive in the original order or $\sqrt{\mathrm{Number \: of \: Walsh \: modes}}$ indices apart (i.e. a supposed order along the columns of the indexes matrix), due to symmetries in their profile.
The random permutation was obtained by using the function ``randperm'' in MatLab.

The parameters of the simulation were: light wavelength of 488 nm, a core diameter of 50 $\mathrm{\mu m}$, a core refractive index of 1.45, NA of 0.22 and a fibre length of 1.5 cm. The window size was set to 62.5 $\mathrm{\mu m}$ x 62.5 $\mathrm{\mu m}$ (2.5 x core radius).

The target of the optimisation was chosen at close to the edge of the fibre core (0.9 x core radius). This was chosen such that the focus location is generated by a larger number of fibre modes.

The SLM pattern for the generation of a focus at the desired location through a straight fibre (inverse radius of 0 $\mathrm{m}^{-1}$) was used as the starting point of the optimisation process. The correction was repeated for a set of inverse bending radii: \{0.5, 1, 2, 3, 4, 5, 7.8, 8.6, 9.2, 10, 12.5, 15, 17.5, 20, 25, 30, 35, 40, 45, 50\} $\mathrm{m}^{-1}$. More samples were used for smaller bends.

Working with the inverse bending radius is easier than with the bending radius itself, as the straight fibre case is represented at 0 $\mathrm{m}^{-1}$ rather than $\infty$ m, thus, facilitating both the calculations involved in the simulation as well as plotting the data.

%subsubs
\subsubsection{Results}

The first step in the simulation process was to observe the achievable intensity if a classic phase-only calibration is used, as described in section \ref{ssec:fibre_calibration}, assuming the amplitude at the SLM is uniform. The SLM is assumed to be in the Fourier plane of the proximal facet. Thus, the propagation from the SLM plane to the proximal facet plane is simulated by using an inverse 2D Fast Fourier Transform (FFT).

First, the TMs for the different bending inverse radii were computed. Starting with a delta function at the distal facet (1 at the focus location and 0 everywhere else) and back propagating through the fibre, using the TM inverse, the desired optical field at the proximal facet was computed. A 2D FFT was then applied to the desired proximal field to find the desired SLM plane optical field. As mentioned previously, the amplitude of every pixel in the SLM plane was set to 1 and the phase information was kept. This is the full correction SLM field.

The first result to be computed was the full correction output intensity of at the target location. This was done by forward propagating the full correction SLM field. 

The first results are in figure \ref{fig:walsh_focus_comparison}

The results are in figure \ref{fig:walsh_full_corr_int}. The intensity varies at different inverse bending radii, however, the variation is on the order of $1\%$, which is most likely due to variations in the amount of light coupling into the fibre and numerical error. These values were later used as normalisation for the optimisation process results.

\begin{sidewaysfigure}
    \centering
    \includegraphics[width=\textwidth]{images/walsh_full_correction_intensity.png}
    \caption{The output intensity of at the target location normalised to the straight fibre case for different inverse bending radii.}
    \label{fig:walsh_full_corr_int}
\end{sidewaysfigure}

Next, the optimisation process for the $64\times64$ size Walsh modes was simulated. This process was repeated for all inverse bending radii. Three examples are given in figure \ref{fig:walsh_opt_4096_individual}. The first apparent observation is that the different permutation perform differently, as expected. furthermore, the general evolution for all permutations is that there is larger improvement at the start of the optimisation process, which reduces over in the latter stages.  

\begin{sidewaysfigure}
    \centering
    \begin{subfigure}[t]{\textwidth}
        \includegraphics[width = \textwidth]{images/walsh_4096_result_invr_2.png}
        \caption{Inverse radius: 2 $\mathrm{m^{-1}}$}
        \label{sfig:walsh_4096_2}
    \end{subfigure}
\end{sidewaysfigure}
\begin{sidewaysfigure}
    \ContinuedFloat
    \begin{subfigure}[t]{\textwidth}
        \includegraphics[width = \textwidth]{images/walsh_4096_result_invr_12p5.png}
        \caption{Inverse radius: 12.5 $\mathrm{m^{-1}}$}
        \label{sfig:walsh_4096_5}
    \end{subfigure}
\end{sidewaysfigure}
\begin{sidewaysfigure}
    \ContinuedFloat
    \begin{subfigure}[t]{\textwidth}
        \includegraphics[width = \textwidth]{images/walsh_4096_result_invr_30.png}
        \caption{Inverse radius: 30 $\mathrm{m^{-1}}$}
        \label{sfig:walsh_4096_12p5}
    \end{subfigure}   
    \caption{Intensity evolution at the target output location normalised to the full correction value, over the Walsh mode optimisation process's iterations. Results for different inverse radii.}
    \label{fig:walsh_opt_4096_individual}
\end{sidewaysfigure}

Another value of interest which can be plotted, is the enhancement factor. This value is defined as the ratio between the output pixel intensity and the total energy, being a good measure of the focus quality. Figure \ref{sfig:walsh_enh_norm} plots the resulting enhancement factor and normalised pixel intensity as a function of iteration number. It is clear that the plot approximately follows the same envelope, but with a variable scaling. Thus, in figure \ref{sfig:walsh_enh_norm_ratio} the ratio between the enhancement factor and the normalised pixel intensity is shown. In terms of what this ratio means, it is $\frac{\mathrm{pixel \ intensity}}{\mathrm{total \ energy}} / \frac{\mathrm{pixel \ intensity}}{\mathrm{fully \ corrected \ pixel \ intensity }}$. Thus, as the fully corrected pixel intensity is a constant and the variable pixel intensity terms cancel, the ratio represents an inversely-proportional envelope of the total energy in the output optical field. As this value increases with iteration number, is signifies that with added Walsh mode modulation, the total energy in the output field decreases up to a point where it plateaus. Thus, the amount of light which couples into the fibre decreases with added modulation. This is an important results to consider, as, while the focus quality increases with respect to the background, the total amount of light coupled into the fibre decreases.

\begin{sidewaysfigure}
    \begin{subfigure}[t]{\textwidth}
        \includegraphics[width = \textwidth]{images/walsh_enhancement_vs_norm_pixel.png}
        \caption{Plots of enhancement ratio and the normalised output pixel intensity.}
        \label{sfig:walsh_enh_norm}
    \end{subfigure}
\end{sidewaysfigure}
\begin{sidewaysfigure} \label{fig:walsh_enh_vs_norm}
   \ContinuedFloat
    \begin{subfigure}[t]{\textwidth}
        \includegraphics[width = \textwidth]{images/walsh_enhancement_vs_norm_pixel_ratio.png}
        \caption{Plots of the ratio between the enhancement ratio and the normalised output pixel intensity.}
        \label{sfig:walsh_enh_norm_ratio}
    \end{subfigure}
    \caption{Comparison between the enhancement ratio and the normalised output pixel intensity. Data generated for a inverse bending radius of 18.25 $m^{-1}$ using an optimisation with 4096 iterations of 64x64 size Walsh modes. First subfigure plots both evolutions as a function of iteration number. The second subplot, plots the ratio between the enhancement factor and iteration number as a function of iteration number.}
\end{sidewaysfigure}


Figure \ref{fig:walsh_opt_final_int} shows the final intensity achieved at the output target location after the optimisation process. To better illustrate the amount of correction, the initial intensity was also plotted with a dashed black line. As expected, the initial non-corrected intensity decreased with a larger bend. However, the optimisation procedure had better results for larger bends. In the worst case investigated the normalised intensity increased from almost $0\%$ to $20\%$ and $45\%$ for the 4096 and 16384 iterations cases, respectively. This is a drastic improvement. Furthermore, it seems that the amount of improvement plateaued, thus, larger and larger bends do not negatively impact the final result of the optimisation procedure.

Another observation is that while continuing with the optimisation process by re-iterating through the permutation did increase the final intensity result, the change is not as much as using larger Walsh modes. This is to be expected, as larger Walsh modes have higher dimensionality, which covers the SLM field space better. Thus, if a larger number of iterations are necessary, it is better to increase the size of the Walsh modes, rather than repeat the iterations with smaller Walsh modes. Though, this is limited by the size of the SLM. 

A final observation is that the original permutation performed worse, while the other two permutations performed similarly. Thus, while the diagonal permutation might not be the best option, and an investigation into which is the best permutation might be desirable, it is a good choice.

\begin{sidewaysfigure}
    \centering
    \begin{subfigure}[t]{\textwidth}
        \includegraphics[width = \textwidth]{images/walsh_4096_final_intensity_all.png}
        \caption{4096 iterations of the 64x64 Walsh modes}
        \label{sfig:walsh_4096_final_64_64_4096}
    \end{subfigure}
\end{sidewaysfigure}
\begin{sidewaysfigure}
    \ContinuedFloat
    \begin{subfigure}[t]{\textwidth}
        \includegraphics[width = \textwidth]{images/walsh_4096_final_intensity_all_4_repeats.png}
        \caption{16384 iterations of the 64x64 Walsh modes}
        \label{sfig:walsh_4096_final_64_64_16384}
    \end{subfigure}
\end{sidewaysfigure}
\begin{sidewaysfigure}
    \ContinuedFloat
    \begin{subfigure}[t]{\textwidth}
        \includegraphics[width = \textwidth]{images/walsh_16384_final_intensity_all.png}
        \caption{16384 iterations of the 128x128 Walsh modes}
        \label{sfig:walsh_16384_final_128_128_16384}
    \end{subfigure}
    \caption{Final output intensity after the optimisation normalised to full correction intensity as a function of inverse bending radius. Dashed black line represents the intensity at the start of the operation.}
    \label{fig:walsh_opt_final_int}
\end{sidewaysfigure}

Finally, a comparison of the foci before bending, after bending with no correction and after correction is given in figure \ref{fig:walsh_focus_comparison} for the case of optimising with Walsh modes of size 128x128 over 16384 iterations at an inverse radius of 18.25 $m^{-1}$.

\begin{sidewaysfigure}
    \centering
    \includegraphics[width=\textwidth]{images/walsh_focus_comparison.png}
    \caption{The achieved focus with a straight fibre, with a bent fibre but no correction (same input as in the straight fibre case) and after 16384 iterations of the Walsh correction algorithm with a Walsh mode size of 128x128. The intensity is in arbitrary units and only a crop around the focus is shown.}
    \label{fig:walsh_focus_comparison}
\end{sidewaysfigure}
%subs

\subsection{Experimental results}

The Walsh mode optimisation process was implemented into the optical system described in section \ref{sssec:optical_system}.

The modulation on the SLM is being controlled thorough the use of a graphical processing unit (GPU), through the use of a Cuda and OpenGL program. Thus, loading the entire Hadamard matrix into the memory of the GPU would be inefficient, as this matrix increases quadratically with respect to the size of the Walsh modes. Thus, each Walsh mode was generated as needed using a Cuda kernel, taking full advantage of the parallelisation strength of a GPU.

Starting from the recursive process of generating a Hadamard matrix, a specific formula for an element of the Hadamard matrix of size $2^m$x$2^m$ is found \cite{sao_hadamard_wiki} in equation \ref{eq:hadamard_element}, where $i$ and $j$ are the row and column indices of the element, $i_k$ and $j_k$ are the bits in the bit-wise (base 2) representation of $i$ and $j$ and $k$ is a bit order index.

\begin{equation} \label{eq:hadamard_element}
    (\mathbf{H}_{2^m})_{i, j} = (-1)^{\sum_k i_k j_k}
\end{equation}

For example, for element (7,12) [element (0111, 1100) in base 2] in a 16x16 Hadamard matrix, $\sum_k i_k j_k = 0\cdot1 + 1\cdot1 + 1\cdot0 + 1\cdot0 = 1$, thus, the element is $(-1)^1 = -1$.

In terms of the algorithm, knowing the index of the Walsh mode generated (j) and the pixel index (i), which can be found from the 2D coordinates of the pixel, the value its associated element in the Walsh mode can be readily computed. The kernel was implemented to compute the value of each Walsh mode element in parallel.

The experiment started by calibrating the fibre and choosing a target pixel. The SLM field which generates the focus, calculated based on the calibration results, was used as a base phase modulation field to which the Walsh mode phase modulation was added. The calibration was done to a subpar standard in order to start from an unoptimised focus. The experiment only used the $64 \times 64$ size Walsh modes.

The first experiment was to check the sinusoidal evolution of a target output intensity as a function of a single Walsh mode coefficient. The Walsh mode chosen had index 5, however, this evolution appears for any Walsh mode, except the first, which represents only a global change in phase and leads to no change in amplitude (technically, this is a sinusoidal variation with amplitude 0).The results are shown in figure \ref{fig:walsh_exp_sinusoidal}. While there is some noise, the sinusoidal variation is apparent.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{images/LV_single_Walsh_single_pixel_intensity_good_plot.png}
    \caption{The output intensity of at the target location as a function of a single Walsh mode coefficient. The chosen Walsh mode had index 5.}
    \label{fig:walsh_exp_sinusoidal}
\end{figure}

The second experiment was to apply the Walsh mode optimisation procedure. The modes were iterated in the original order. At every iteration, in order to limit noise, three acquisitions are made for each measurement and the result used is the average value. Unfortunately, due to instabilities in the optical system control software, only 100 iterations were able to be computed. However, even 100 iterations clearly show an improvement in intensity, whose rate decreases with larger number of iterations. An observation to be made is that in the real system, unlike the simulations, the intensity does not increase at every single iteration. This was most likely due to noise, thus more measures would be needed to minimise it, especially as the number of iterations increases, as the amount of improvement decreases. Though, at some point the grayscale resolution of the camera will become a limitation.

\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{images/LV_Walsh_optimisation_4096_modes_4096_iter_100_good_plot.png}
    \caption{The evolution of output intensity of at the target location over the iterative Walsh optimisation process.}
    \label{fig:walsh_exp_optimisation}
\end{figure}

% subs
\subsection{Conclusion}

The sensorless adaptive optics approach provides a means to recovering focusing capabilities at a desired location at or near the distal facet plane. The approach only requires that the intensity of the focus at the target be measured and the input field be modulated in phase in its Fourier plane through the use of a modulator such as an SLM. The approach was shown, through simulations, to be able to recover a large factor of the lost intensity in the focus due to bending of the optical fibre. However, in order to do so, a large number of iterations are needed, which depending on the application might not be desirable, as the time required for the optimisation to be completed is long.

In conclusion, while it has some limitations, the sensorless approach to recovering a focus at the output of the fibre is a feasible approach, though, it could be improved or replaced with a faster method. 

\end{document}