Phase retrieval

Schematic view of the error reduction algorithm for phase retrieval

The error reduction is a generalization of the Gerchberg–Saxton algorithm. It solves for ${\displaystyle f(x)}$ from measurements of ${\displaystyle |F(u)|}$ by iterating a four-step process. For the ${\displaystyle k}$th iteration the steps are as follows:

Step (1): ${\displaystyle G_{k}(u)}$, ${\displaystyle \phi _{k}}$, and ${\displaystyle g_{k}(x)}$ are estimates of, respectively, ${\displaystyle F(u)}$, ${\displaystyle \psi }$ and ${\displaystyle f(x)}$. In the first step we calculate the Fourier transform of ${\displaystyle g_{k}(x)}$:

${\displaystyle G_{k}(u)=|G_{k}(u)|e^{i\phi _{k}(u)}={\mathcal {F}}(g_{k}(x)).}$

Step (2): The experimental value of ${\displaystyle |F(u)|}$, calculated from the diffraction pattern via the signal equation, is then substituted for ${\displaystyle |G_{k}(u)|}$, giving an estimate of the Fourier transform:

${\displaystyle G'_{k}(u)=|F(u)|e^{i\phi _{k}(u)},}$

where the ' denotes an intermediate result that will be discarded later on.

Step (3): the estimate of the Fourier transform ${\displaystyle G'_{k}(u)}$ is then inverse Fourier transformed:

${\displaystyle g'_{k}(x)=|g'_{k}(x)|e^{i\theta '_{k}(x)}={\mathcal {F}}^{-1}(G'_{k}(u)).}$

Step (4): ${\displaystyle g'_{k}(x)}$ then must be changed so that the new estimate of the object, ${\displaystyle g_{k+1}(x)}$, satisfies the object constraints. ${\displaystyle g_{k+1}(x)}$ is therefore defined piecewise as:

${\displaystyle g_{k+1}(x)\equiv {\begin{cases}g'_{k}(x),&x\notin \gamma ,\\0,&x\in \gamma ,\end{cases}}}$

where ${\displaystyle \gamma }$ is the domain in which ${\displaystyle g'_{k}(x)}$ does not satisfy the object constraints. A new estimate ${\displaystyle g_{k+1}(x)}$ is obtained and the four step process is repeated.

This process is continued until both the Fourier constraint and object constraint are satisfied. Theoretically, the process will always lead to a convergence,[1] but the large number of iterations needed to produce a satisfactory image (generally >2000) results in the error-reduction algorithm by itself being unsuitable for practical applications.

The hybrid input-output algorithm is a modification of the error-reduction algorithm - the first three stages are identical. However, ${\displaystyle g_{k}(x)}$ no longer acts as an estimate of ${\displaystyle f(x)}$, but the input function corresponding to the output function ${\displaystyle g'_{k}(x)}$, which is an estimate of ${\displaystyle f(x)}$.[1] In the fourth step, when the function ${\displaystyle g'_{k}(x)}$ violates the object constraints, the value of ${\displaystyle g_{k+1}(x)}$ is forced towards zero, but optimally not to zero. The chief advantage of the hybrid input-output algorithm is that the function ${\displaystyle g_{k}(x)}$ contains feedback information concerning previous iterations, reducing the probability of stagnation. It has been shown that the hybrid input-output algorithm converges to a solution significantly faster than the error reduction algorithm. Its convergence rate can be further improved through step size optimization algorithms.[2]

${\displaystyle g_{k+1}(x)\equiv {\begin{cases}g'_{k}(x),&x\notin \gamma ,\\g_{k}(x)-\beta {g'_{k}(x)},&x\in \gamma .\end{cases}}}$

Here ${\displaystyle \beta }$ is a feedback parameter which can take a value between 0 and 1. For most applications, ${\displaystyle \beta \approx 0.9}$ gives optimal results.{Scientific Reports volume 8, Article number: 6436 (2018)}

For a two dimensional phase retrieval problem, there is a degeneracy of solutions as ${\displaystyle f(x)}$ and its conjugate ${\displaystyle f^{*}(-x)}$ have the same Fourier modulus. This leads to "image twinning" in which the phase retrieval algorithm stagnates producing an image with features of both the object and its conjugate.[3] The shrinkwrap technique periodically updates the estimate of the support by low-pass filtering the current estimate of the object amplitude (by convolution with a Gaussian) and applying a threshold, leading to a reduction in the image ambiguity.[4]

Given a discrete signal ${\displaystyle \mathbf {x} =(f[0],f[1],...,f[N-1])^{T}}$ which is sampled from ${\displaystyle f(x)}$. We use a window of length W: ${\displaystyle \mathbf {w} =(w[0],w[1],...,w[W-1])^{T}}$ to compute the STFT of ${\displaystyle \mathrm {f} }$, denoted by ${\displaystyle \mathbf {Y} }$:

${\displaystyle Y[m,r]=\sum _{n=0}^{N-1}{f[n]w[rL-n]e^{-i2\pi {\frac {mn}{N}}}}}$

for ${\displaystyle 0\leq m\leq N-1}$ and ${\displaystyle 0\leq r\leq R-1}$, where the parameter ${\displaystyle L}$ denotes the separation in time between adjacent short-time sections and the parameter ${\displaystyle R=\left\lceil {\frac {N+W-1}{L}}\right\rceil }$ denotes the number of short-time sections considered.

The other interpretation (called sliding window interpretation) of STFT can be used with the help of discrete Fourier transform (DFT). Let ${\displaystyle w_{r}[n]=w[rL-n]}$ denotes the window element obtained from shifted and flipped window ${\displaystyle \mathbf {w} }$. Then we have

${\displaystyle \mathbf {Y} =[\mathbf {Y} _{0},\mathbf {Y} _{1},...,\mathbf {Y} _{R-1}]}$, where ${\displaystyle \mathbf {Y} _{r}=\mathbf {x} \circ \mathbf {w} _{r}}$.

Let ${\displaystyle {Z}_{w}[m,r]=|Y[m,r]|^{2}}$ be the ${\displaystyle N\times R}$ measurements corresponding to the magnitude-square of the STFT of ${\displaystyle \mathbf {x} }$, ${\displaystyle \mathbf {W} _{r}}$ be the ${\displaystyle N\times N}$ diagonal matrix with diagonal elements ${\displaystyle \left(w_{r}[0],w_{r}[1],\ldots ,w_{r}[N-1]\right).}$ STFT phase retrieval can be stated as:

Find ${\displaystyle \mathbf {x} }$ such that ${\displaystyle Z_{w}[m,r]=\left|\left\langle \mathbf {f} _{m},\mathbf {W} _{r}\mathbf {x} \right\rangle \right|^{2}}$ for ${\displaystyle 0\leq m\leq N-1}$ and ${\displaystyle 0\leq r\leq R-1}$, where ${\displaystyle \mathbf {f} _{m}}$ is the ${\displaystyle m}$-th column of the ${\displaystyle N}$-point inverse DFT matrix.

Intuitively, the computational complexity growing with ${\displaystyle N}$ makes the method impractical. In fact, however, for the most cases in practical we only need to consider the measurements corresponding to ${\displaystyle 0\leq m\leq M}$, for any parameter ${\displaystyle M}$ satisfying ${\displaystyle 2W\leq M\leq N}$.

To be more specifically, if both the signal and the window are not vanishing, that is, ${\displaystyle x[n]\neq 0}$ for all ${\displaystyle 0\leq n\leq N-1}$ and ${\displaystyle w[n]\neq 0}$ for all ${\displaystyle 0\leq }$ ${\displaystyle n\leq W-1}$, signal ${\displaystyle \mathbf {x} }$ can be uniquely identified from its STFT magnitude if the following requirements are satisfied:

1. ${\displaystyle L<W\leq N/2}$,
2. ${\displaystyle 2W\leq M\leq N}$.

The proof can be found in Jaganathan' s work,[5] which reformulates STFT phase retrieval as the following least-squares problem:

${\displaystyle \min _{\mathbf {x} }\sum _{r=0}^{R-1}\sum _{m=0}^{N-1}\left(Z_{w}[m,r]-\left|\left\langle \mathbf {f} _{m},\mathbf {W} _{r}\mathbf {x} \right\rangle \right|^{2}\right)^{2}}$.

The algorithm, although without theoretical recovery guarantees, empirically able to converge to the global minimum when there is substantial overlap between adjacent short-time sections.

To establish recovery guarantees, one way is to formulate the problems as a semidefinite program (SDP), by embedding the problem in a higher dimensional space using the transformation ${\displaystyle \mathbf {X} =\mathbf {x} \mathbf {x} ^{\ast }}$ and relax the rank-one constraint to obtain a convex program. The problem reformulated is stated below:

Obtain ${\displaystyle \mathbf {\hat {X}} }$ by solving:

$${\displaystyle {\begin{aligned}&\mathrm {minimize} ~~~\mathrm {trace} (\mathbf {X} )\\[6pt]&\mathrm {subject~to} ~~Z[m,r]=\mathrm {trace} (\mathbf {W} _{r}^{\ast }\mathbf {f} _{m}\mathbf {f} _{m}^{\ast }\mathbf {W} _{r}\mathbf {X} )\\[0pt]&~~~~~~~~~~~~~~~~~~~\mathbf {X} \succeq 0\end{aligned}}}$$

for ${\displaystyle 1\leq m\leq M}$ and ${\displaystyle 0\leq r\leq R-1}$

Once ${\displaystyle \mathbf {\hat {X}} }$ is found, we can recover signal ${\displaystyle \mathbf {x} }$ by best rank-one approximation.