## Abstract

The gyrator transform is a useful tool for optical information processing applications. In this work we propose a multi-stage phase retrieval approach based on this operation as well as on the well-known Gerchberg-Saxton algorithm. It results in an iterative algorithm able to retrieve the phase information using several measurements of the gyrator transform power spectrum. The viability and performance of the proposed algorithm is demonstrated by means of several numerical simulations and experimental results.

©2009 Optical Society of America

## 1. Introduction

Phase retrieval plays an important role in science, engineering, and in particular in optical imaging and information processing. In general, its experimental application demands a flexible and robust setup which is an important issue. There exist several approaches for phase retrieval used for example in the estimation of the depth dimension of objects, characterization of optical elements, etc. Since holography permits the recording of the complex field amplitude, some of these phase retrieval approaches exploit interferometric techniques [1, 2, 3, 4], which requires the use of a reference signal. Alternatively, phase retrieval approaches based only on the measurement of intensity distributions during signal propagation are applied.

These methods are either iterative or non-iterative involving digital post processing. In particular, the non-iterative approaches are able to retrieve the phase information by means of transforming Guigay equations [5, 6], by solving the so-called transport of intensity equation [7], etc. Nevertheless, they are often mathematically complex. Other option is to use iterative techniques based on the well-known Gerchberg-Saxton (GS) algorithm [8]. The GS algorithm is widely used in different optical applications and permits to retrieve the input phase from the known amplitude distributions corresponding to input and output planes of the optical setup. Its main problem for phase retrieval is the phase ambiguity resulting with converging to local minima solutions, which yields an incorrect phase distribution. There exist several works trying to solve this problem by using statistical methods [9], but they demand a high computational complexity. Blind de-convolution methods [10] speed up the convergence rate but do not solve the local minima problem either.

Other iterative approaches have been also proposed, for instance the gradient search method [11]. However, unlike the GS approach, such techniques are only suitable for special conditions. Alternatively, there exist methods in microscopy using plurality of de-focused images [12, 13]. These methods use de-focused images as well as a focused image to reconstruct the phase distribution of one of the planes, by using the known relations between the images. Later on, other works were introduced combining this idea with an iterative GS algorithm [14, 15, 16]. This approach was found to be faster than the iterative GS technique and it also solved the phase ambiguity issue due to the usage of multiple constraint planes. Originally GS technique involved a Fourier or Fresnel transform relation between the input and output planes. Nevertheless, later on other transformations relating such planes were suggested [17] as well.

In this work we apply a GS-based approach involving several intensity distributions working as constraint planes. In particular, we consider the gyrator transform (GT) which is mathematically defined as a linear canonical integral transform that produces the twisted rotation in position-spatial frequency planes of phase space [18]. Thus, the GT operation of a two-dimensional function *f _{i}*(

**r**

_{i}) (complex field amplitude) with parameter

*α*, known as transformation angle, is given by

where **r**
* _{i,o}* = (

*x*) are the input and output coordinates, respectively. This transformation can be implemented optically by using cylindrical lenses which are properly rotated to change the transformation angle

_{i,o}, y_{i,o}*α*[19]. The experimental implementation of the GT has been demonstrated in Ref. [18].

In our approach at each iteration loop an iterative calculation based on the GT operation for a certain set of angles *α* is performed. After calculating all the considered angles, an average estimation is obtained yielding a retrieved phase to be used as input for the next iteration loop. The validity of the method is demonstrated by numerical simulations as well as by experimental results.

We remark that GT operation for different angles *α* can be performed by only the proper rotation of cylindrical lenses, where the distance between them and input-output planes are fixed. This fact avoid several problems arising from misalignment which is an important issue in optical systems that require axial movements. Therefore the proposed approach results attractive for experimental phase retrieval applications.

## 2. Multi-stage phase retrieval algorithm

The proposed phase retrieval algorithm is based on the measurements of the intensity distributions associated with GT power spectrum |*F ^{α}*(

**r**

_{o})|

^{2}for various angles

*α*and |

*f*(

_{i}**r**

*)|*

_{i}^{2}= |

*F*

^{0}(

**r**

_{i})|

^{2}. Note that |

*F*(

^{π}**r**

_{o})|

^{2}= |

*F*(-

_{i}**r**

_{o})|

^{2}. Furthermore we will use the following notations for the input signal’s amplitude |

*F*(

_{i}**r**

_{i})| =

*E*(

_{in}*x*,

*y*) and the amplitude of its GT operation |

*F*(

^{αj}**r**

_{o})| =

*E*(

_{out}^{j}*u*,

*v*), where

*j*= 1, 2, …,

*M*.

Let us now explain the flow chart of the proposed approach as it is depicted in Fig. 1. An initial random phase *ϕ*
_{1}(*x*,*y*) ∈ (0, 2*π*) is chosen. The retrieved phase distribution *ϕ _{k}*(

*x*,

*y*) at iteration number

*k*is obtained as follows:

- The phase estimation from the previous iteration step is multiplied by the input amplitude distribution
*E*(_{in}*x*,*y*), resulting with*E*(_{in}*x*,*y*) exp [*iϕ*(_{k}*x*,*y*)]. - Perform the following calculation
*M*times:- Gyrator transform of
*E*(_{in}*x*,*y*) exp [*iϕ*(_{k}*x*,*y*)] with the transformation angle*α*resulting in a function denoted by:_{j}*B*(_{k}^{j}*u*,*v*) exp*E*(_{out}^{j}*u*,*v*) . - Replace amplitude
*B*(_{k}^{j}*u*,*v*) with*E*(_{out}^{j}*u*,*v*). - Gyrator transform of
*E*(_{out}^{j}*u*,*v*) exp [*iϕ*(_{k}^{j}*u*,*v*)] with the transformation angle -*α*(which equals to the inverse GT) resulting with_{j}*A*(_{k}^{j}*x*,*y*) exp [*iϕ*(_{k}^{j}*x*,*y*)]. - Use the retrieved phase
*ϕ*(_{k}^{j}*x*,*y*) for the next iteration corresponding to angle*α*_{j+1}.

The previous calculation leads to a new phase estimation to be used at the next loop: *ϕ*
_{k+1}(*x*,*y*) = *ϕ _{k}^{M}*(

*x*,

*y*). The phase retrieval algorithm ends after

*N*loops yielding a retrieved phase

*ϕ*(

_{N}*x*,

*y*). Therefore, this approach requires

*N*×

*M*iterations, where

*M*is the number of transformation angles

*α*or considered constrain planes.

In order to demonstrate the phase retrieval capabilities of the proposed approach, we preformed numerical simulations following the parameters used further in the experimental GT setup. In all simulations we have chosen the parameters to match the values of the experimental realization for the GT, which is described by:

where *λ* is the wavelength, *z* corresponds to fixed axial interval between the lenses and *s*
^{2} = 2*λz* is a constant coordinate normalization [18]. In our experimental setup *λ* = 532 nm and *z* = 0.94 cm.

As input signals we consider the Hermite-Gaussian (HG) modes

and Laguerre-Gaussian (LG) modes

where: *r*
^{2} = *x*
^{2} + *y*
^{2}, *H _{m}* is the Hermite polynomial, $w=\sqrt{2\mathrm{\lambda z}}$ is the beam waist, and

*L*is the Laguerre polynomial with radial index

_{p}^{l}*p*and azimuthal index

*l*. In particular we study the HG

_{4,3}and LG

_{3,1}

^{+}mode, see Fig. 2. We underline that HG

_{4,3}is transformed into LG

_{3,1}

^{+}mode under the GT for

*α*= (2

*q*+ 1)

*π*/4 and vice versa, where

*q*is an integer,

*p*= min(

*m*,

*n*) and

*l*= |

*m*-

*n*|.

Now we aim to demonstrate that using several (rather than a single) constraint planes for different transformation angles of the GT, produce improved phase retrieval. In particular we study the phase-retrieval performance corresponding to HG_{4,3} mode by using one, two and three constraint planes (*M* = 1, 2, 3…) in which the amplitude distribution is known. The retrieved phase distributions after 30 iterations for HG_{4,3} mode are displayed in Fig. 3(a) with constraint planes involving angles *α*
_{1} = 45^{o}, *α*
_{2} = 65^{o} and *α*
_{3} = 80^{o}. Notice that the GS algorithm allows to reconstruct the phase up to a constant. Because of that, the phase distributions reconstructed from two and three constraint planes differ by a constant phase, see Fig. 2(a).

The relative root mean square (RMS) error, also known as Abserr, of the retrieved phase at each iteration is given by

where *ϕ*(*x*,*y*) is the phase to retrieve and *ϕ _{rt}*(

*x*,

*y*) is the estimated phase distribution at each iteration.

The evolution of RMS error for the case of HG_{4,3} mode is displayed in Fig. 3(b), respectively calculated for one, two and three constraint planes. We conclude that increasing the number of constraint planes from one to two and three, speeds up the converging of the RMS error as one appreciates in Fig. 3(b). This fact has also been proven for other values of the transformation angles yielding a retrieved phase with a RMS error value of 0.01 or less.

In order to demonstrate the RMS error evolution as a function of the number of constraint planes with low number of iterations, we now consider LG^{+}
_{3,1} mode as input signal. Specifically, we use only 36 iterations and constraint planes corresponding to the GT operation at angles *α* = 45^{o}, 55^{o}, 65^{o}, 70^{o} and 80^{o}. The results displayed in Fig. 4, demonstrate that the phase distribution is successfully retrieved by using at least three constraint planes obtaining a RMS error value of 0.006 after 36 iterations. In conclusion, when a low number of iterations is utilized more constraint planes are needed in order to obtain an accurate retrieved phase (low RMS value). Meanwhile, as it was previously demonstrated, the RMS error can be reduced to 0.01 or less for the case of only two constraint planes by increasing the number of iterations.

To simulate the phase retrieval from measurements of noisy intensity distributions, we consider a normal random additive noise. In particular, we produce this random noise yielding a signal-to-noise ratio (SNR) of 20 dB which is added to each power spectrum or constraint plane. We have chosen this noise and SNR value because it is similar to the one observed in our experimental data. For instance, the constraint plane *α*= 80^{o} displayed in Fig. 5(a) is now disturbed by this additive noise as it is shown in Fig. 5(b). Here we perform the same calculation as in the noiseless case: five constraint planes corresponding to angles *α* = 45^{o}, 55^{o}, 65^{o}, 70^{o} and 80^{o}. As it is demonstrated in Fig. 6 the phase distribution is successfully retrieved even in this noisy case by using at least three constraint planes and after 36 iterations (RMS of 0.02). Therefore the proposed algorithm is comparatively robust for noise, SNR of 20 dB in our case. We have also studied the case of high-moderate noisy images. For example in the Table 1 the RMS error values obtained in the case of SNR of 15 dB are also displayed. In contrast to the noiseless case (see Table 1) now more constraint planes are needed in order to obtain a similar RMS error value, as it was expected. We underline that the RMS error is again minimized by increasing the number of iterations as well as constraint planes utilized until the convergence occurs.

## 3. Experimental results

The GT setup is constructed by using three generalized lenses (with L_{3} = L_{1}) and two fixed intervals *z*, see Fig. 7(a). Each generalized lens is an assembled set of two thin convergent cylindrical lenses, Fig. 7(b). The rotation angles of the cylindrical lenses are *ϕ*
_{1} = -*φ*
_{1,2} and *ϕ*
_{2} = - (*ϕ*
_{1} + *π*/2), where the angles *φ*
_{1} and *φ*
_{2} correspond to generalized lens L_{1} and L_{2} as follows: sin2*φ*
_{1} = cot (*α*/2) /2 and sin2*φ*
_{2} = sin *α*. Notice that we have chosen a coordinate normalization of *s*
^{2} = 2*λz* instead of the one reported in Ref. [18] (*s*
^{2} = *λz*).

Since this application only demands the acquisition of the intensity distribution then the third generalized lens is not required. Therefore the experimental setup consists of two generalized lenses (four cylindrical lenses), see Fig. 8(a). We remark that each generalized lens is an assembled set of two identical plano-convex cylindrical lenses which are fabricated from N-BK7 glass. Moreover, in order to obtain an accurate assembly the cylindrical lenses were set into a cage-system with rotation mount, see Fig. 8(b).

Let us consider the LG^{+}
_{3,1} mode [see Fig. 9(a)] as input function in order to demonstrate experimentally the phase-retrieval performance. Intensity and phase distributions of its GT for *α* = 135^{o}, 115^{o} and 100^{o} are displayed in Fig. 9, first and second rows, respectively. Here we perform the phase retrieval algorithm to the experimental intensity distributions, shown at the third row in Fig. 9. The phase distribution is successfully retrieved with a RMS error of 0.30 after 45 iterations by using three constrain planes, see Fig. 10.

We underline that there exists a wavefront distortion arising from the experimental setup, as it is observed in the intensity distributions displayed at the third row in Fig. 9. Therefore the retrieved phase shown in Fig. 10 differs from the theoretical one, see Fig. 2(b), by an additional phase distribution *W*(*x*,*y*) associated to the wavefront distortion. This additional phase signal can be extracted from the retrieved phase *R*(*x*,*y*) by using the phase distribution of the ideal LG^{+}
_{3,1} mode considered as input signal: *W*(*x*,*y*) = *R*(*x*,*y*) - arg [LG^{+}
_{3,1}], see Fig. 11.

It is expected that in the case of exact knowledge of this additional phase *W*(*x*,*y*) the residual RMS errors have to be smaller than it is presented in Fig. 10, where only the ideal phase of LG^{+}
_{3,1} was used as reference for the RMS convergence analysis. Besides, the RMS convergence is expected to be similar to the one previously studied for the ideal case, Fig. 4. In contrast with the simulations corresponding to a distortion-free case (Fig. 4), now the RMS convergence is reached almost at the same number of iterations with RMS values of 0.35 and 0.30 for the case of two and three constraint planes respectively, see Fig. 10. The RMS evolution in this case (Fig. 10) has to be understood as a measure of the difference between the phase distribution experimentally retrieved and the one corresponding to an ideal case in which the optical setup as well as the input signal are free of such distortions.

Notice that in our case the input signal is generated by using two spatial light modulators (SLM), see Ref [18]. The nonlinear response of the SLM and the deviation from flatness of its reflective surface lead to distortions in phase modulation that also contributes to the overall wavefront distortion. We remind that the experimental data are also affected by an additive random noise arising from the CCD camera, which does not imply significant distortions in the retrieved phase. This effect was studied and demonstrated in the previous section by means of numerical simulation for SNR of 20 dB, see Fig. 6. Nevertheless, the experimental results are in good agreement with the theoretical predictions demonstrating the algorithm feasibility.

## 4. Conclusions

In this paper we have presented an iterative phase retrieval algorithm based upon the gyrator transform (GT) which applies several, rather than a single, planes with amplitude constrains, along the iterative process. We have demonstrated that increasing the number of constraint planes speeds up the algorithm convergence rate. It also permits to obtain an accurate retrieved phase distribution (RMS of 0.02) even in the case of noisy images with SNR of 20 dB. The algorithm have been applied to the experimental data. The results are in good agreement with the theoretical predictions demonstrating the practical feasibility of the proposed method.

We underline that the GT optical setup proposed for the measurements of the required GT power spectra is a flexible and robust system that does not require axial movements. The variation of the GT parameter is achieved by only the proper rotation of cylindrical lenses. These facts make attractive the application of the discussed approach for phase retrieval of optical beams.

## Acknowledgments

The financial support of the Spanish Ministry of Science and Innovation under project TEC2008-04105 and the Santander-Complutense project PR-34/07-15914 are acknowledged. José A. Rodrigo gratefully thanks a “Juan de la Cierva” grant.

## References

**1. **D. Gabor, “A new microscopic principle,” Nature **161**, 777–778 (1948). [CrossRef] [PubMed]

**2. **E. N. Leith and J. Upatnieks, “Reconstructed wavefronts and communication theory,” J. Opt. Soc. Am. **52**, 1123–1128 (1962). [CrossRef]

**3. **E. N. Leith and J. Upatnieks, “Wavefront reconstruction with continuous-tone objects,” J. Opt. Soc. Am. **53**, 1377–1381 (1963). [CrossRef]

**4. **E. N. Leith and J. Upatnieks, “Wavefront reconstruction with diffused illumination and three-dimensional objects,” J. Opt. Soc. Am. **54**, 1295–1301 (1964). [CrossRef]

**5. **J. P. Guigay, “Fourier-transrorm analysis of Fresnel diffraction patterns and in-line holograms,” Optik **49**, 121–125 (1977).

**6. **M. R. Teague, “Deterministic phase retrieval: a Green’s function solution,” J. Opt. Soc. Am. **73**, 1434–1441 (1983). [CrossRef]

**7. **T. E. Gureyev and K. A. Nugent, “Rapid quantitative phase imaging using the transport of intensity equation,” Opt. Commun. **133**, 339–346 (1997). [CrossRef]

**8. **R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of phase from image and diffraction plane pictures,” Optik **35**, 237–246 (1972).

**9. **M. Nieto-Vesperinas, R. Navarro, and F. J. Fuentes, “Performance of a simulated-annealing algorithm for phase retrieval,” J. Opt. Soc. Am. A **5**, 30–38 (1988). [CrossRef]

**10. **J. H. Seldin and J. R. Fienup, “Iterative blind deconvolution algorithm applied to phase retrieval,” J. Opt. Soc. Am. A **7**, 428–433 (1990). [CrossRef]

**11. **J. R. Fienup, “Phase-retrieval algorithms for a complicated optical system,” Appl. Opt. **32**, 1737–1746 (1993). [CrossRef] [PubMed]

**12. **D. L. Misell, “A method for the solution of the phase problem in electron microscopy,” J. Phys. D: Applied Physics **6**, L6–L9 (1973). [CrossRef]

**13. **D. C. Redding, S. Basinger, A. Lowman, A. Kissil, P. Bely, R. Burg, R. Lyon, G. Mosier, M. Femiano, M. Wilson, R. G. Schunk, L. Craig, D. Jacobson, J. Rakoczy, and J. Hadaway, “Wavefront Sensing and Control for a Next-Generation Space Telescope,” Proc. SPIE **3356**, 758 (1998). [CrossRef]

**14. **D. S. Acton, P. D. Atcheson, M. Cermak, L. K. Kingsbury, F. Shi, and D. C. Redding, “James Webb Space Telescope wavefront sensing and control algorithms,” Proc. SPIE **5487**, 887 (2004). [CrossRef]

**15. **G. R. Brady and J. R. Fienup, “Nonlinear optimization algorithm for retrieving the full complex pupil function,” Opt. Express **14**, 474–486 (2006). [CrossRef] [PubMed]

**16. **B. H. Dean, D. L. Aronstein, J. S. Smith, R. Shiri, and D. S. Acton, “Phase retrieval algorithm for JWST Flight and Testbed Telescope,” Proc. SPIE **6265**, 626511 (2006). [CrossRef]

**17. **Z. Zalevsky, D. Mendlovic, and R. G. Dorsch, “Gerchberg-Saxton algorithm applied in the fractional Fourier or the Fresnel domain,” Opt. Lett. **21**, 842 (1996). [CrossRef] [PubMed]

**18. **J. A. Rodrigo, T. Alieva, and M. L. Calvo, “Experimental implementation of the gyrator transform,” J. Opt. Soc. Am. A **24**, 3135–3139 (2007). [CrossRef]

**19. **J. A. Rodrigo, T. Alieva, and M. L. Calvo, “Optical system design for orthosymplectic transformations in phase space,” J. Opt. Soc. Am. A **23**, 2494–2500 (2006). [CrossRef]