International Journal of Optics

Volume 2018, Article ID 9560401, 10 pages

https://doi.org/10.1155/2018/9560401

## Measuring the Phase of an Optical Field from Two Intensity Measurements: Analysis of a Simple Theoretical Model

Oryx Optics Ltd, 19 The Elms, Athlone, Co. Westmeath, Ireland

Correspondence should be addressed to Damien P. Kelly; moc.liamg@yllekpneimad

Received 21 June 2018; Accepted 16 September 2018; Published 1 October 2018

Guest Editor: Wei Liu

Copyright © 2018 Damien P. Kelly. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Under the scalar paraxial approximation, an optical wavefield is considered to be complex function dependent on position; i.e., at a given location in space the optical field is a complex value with an intensity and phase. The optical wavefield propagates through space and can be modeled using the Fresnel transform. Lenses, apertures, and other optical elements can be used to control and manipulate the wavefield and to perform different types of signal processing operations. Often these optical systems are described theoretically in terms of linear systems theory leading to a commonly used Fourier optics framework. This is the theoretical framework that we will assume in this manuscript. The problem which we consider is how to recover the phase of an optical wavefield over a plane in space. While today it is relatively straightforward to measure the intensity of the optical wavefield over a plane using CMOS or CCD sensors, recovering the phase information is more complicated. Here we specifically examine a variant of the problem of phase retrieval using two intensity measurements. The intensity of the optical wavefield is recorded in both the image plane and the Fourier plane. To make the analysis simpler, we make a series of important theoretical assumptions and describe how in principle the phase information can be recovered. Then, a deterministic but iterative algorithm is derived and we examine the characteristics and properties of this algorithm. Finally, we examine some of the theoretical assumptions we have made and how valid these assumptions are in practice. We then conclude with a brief discussion of the results.

#### 1. Introduction

Visible light is an electromagnetic field with a wavelength that ranges roughly from ultraviolet (300 nm) to infrared (800 nm). An electromagnetic field is vectorial in nature with three electric and magnetic orthogonal field components. When designing a microwave radio receiver (long wave wavelengths relative to visible light) one must consider how these different magnetic and electric field components interact with each other and with conducting strips of metal and dielectric substrates that are used for impedance matching, mixing, and waveguiding. The resulting equations (that follow from Maxwell’s equations) are quite complex and require intensive numerical calculation procedures since the different vectorial components of the EM field interact with each other forming coupled vectorial equations [1].

Fortunately when examining optical problems in the visible regime it is often possible to make the scalar approximation. This means that we can assume that the different vectorial components of the optical wavefield do not “see” each other and hence each vectorial component can be treated independently as a separate scalar problem. This approximation is very useful and quite accurate and can be used to analyze a wide range of important and practical optical problems.

Therefore in this manuscript we assume that a scalar description of light propagation is valid and use the Fresnel transform to relate an optical field in one plane to that in another plane, where the optical planes are separated from each other by an axial distance . The light sources used are assumed to be of a definite monochromatic temporal frequency and are both temporally and spatially coherent [2, 3]. Optical elements used to shape the form of the illuminating wavefield are assumed to be “thin” and operate on an incident optical field, multiplying it by a function describing the optical element over a plane [4]. The assumptions we have described place us in a Fourier optics framework [4]. Fourier optics is an important branch of optical theory since it allows the development of relatively simple and intuitive models of optical systems. These models provide significant insight into the characteristics of the underlying optical systems and are also reasonably accurate. Furthermore Fourier optics builds important bridges between disciplines in particular allowing optical systems to be interpreted in terms of signal processing operations that are commonly used in electrical engineering, communication theory, and control engineering [4–9]. This type of viewpoint can be extended still further with the development of mixed space/time transforms or mixed space/spatial frequency transforms such as the fractional Fourier transform [10], the Linear Canonical Transform (LCT), and perhaps even more generally the form of Wigner distribution functions [11–17].

We note that in this theoretical framework if we know the intensity and phase of the optical field in a give optical plane, it is then possible to calculate (numerically or under certain circumstances analytically) the field distribution in any other plane using the LCT [18, 19] and hence we can image over a 3D volume, with varying levels of spatial resolution [19].

In a practical sense then it is important to be able to measure experimentally both the intensity and the phase of an optical wavefield in a given plane. Nowadays high quality low noise CCD and CMOS sensors are commercially available and widely used. With these devices it is relatively straightforward to measure the intensity of an optical field in a particular plane. We will now briefly discuss several techniques for indirectly measuring the phase distribution of the field.

##### 1.1. Digital Holography

Digital holography (DH) was first proposed by Gabor himself in [20–22] and the first experimental implementation was done by Goodman et al. in [23]; later pioneering work was done by Yaroslawski, [24], Onural and Scott [25], and Schnars and Jüptner [26]. More recent theoretical and experimental analysis of these systems is given [19, 27–36]. In holography we need both an object wave and a so-called reference wave. Light from the laser source in the optical setup is divided by a beam splitter into a reference and illuminating plane waves. The illuminating plane wave is directed by a mirror so that it illuminates the object we wish to examine. The wavefront interacts with a “thin” object and the resulting object wavefield scatters and propagates a distance to the CCD/CMOS plane. Meanwhile the reference wave is directed through the optical system where it is recombined with the object wave at a second beam splitter. Thus at the sensor plane we record the interference between the reference and object wavefields thereby recording a hologram. A drawback of this approach is that the additional optics parts are required; a clean stable reference beam is essential; and several undesired “noise” terms, two DC terms and a twin image term, are produced by the recording process. These terms can be removed by capturing several phase shifted intensity images and performing some postprocessing operations. Single shot off-axis holography with a tilted plane reference wave is possible; however the spatial resolution of the resulting real image term (object wavefield) is significantly reduced. The accuracy and limitations of this approach have already been examined in several publications; see, for example, [33, 37].

##### 1.2. Transport of Intensity

These approaches are quite different in nature and require the capture of two separate intensity distributions. Usually these distributions are captured in an image plane and in a slightly defocused plane. Using the so-called transport of intensity equation, see, for example, [38, 39], it is possible to recover an estimate for the phase distribution. There are several interesting aspects of the reconstruction of the phase; firstly the reconstructed phase is not wrapped over 2 and secondly with some reconstruction approaches, namely, an FFT based approach to indefinite integration, see [40], low frequencies can be more difficult to recover.

##### 1.3. Iterative Phase Retrieval

Finally we turn our attention to iterative phase retrieval approaches. Here usually two or more intensity distributions are captured in different LCT domains, or in physical optics terms the intensity is captured in several different optical planes. Often these different optical planes are separated from each other by sections of free-space and combinations of lenses, but also in the image and corresponding Fourier transform planes. The origin of these types of algorithms dates back to the iterative Gerschberg-Saxton (GS) approach first published in the 70s and was further developed Fienup et al. [41, 42]. The initial algorithms assumed that two measurements were made and the intensity distribution was recorded in the image plane and in the Fourier plane. Later algorithms extended these approaches to work with multiple different Fresnel domains [43–47].

The image-Fourier plane (GS) algorithm works as follows:(1)One starts with an initial guess of the phase distribution in the image plane(2)We then numerically calculate the Fourier transform of this distribution and compare the calculated intensity in the Fourier plane to the experimentally measured result. The calculated intensity distribution will differ from the experimentally measured distribution because our initial guess at the phase distribution will not be correct(3)At this stage we replace the incorrect numerically calculated intensity distribution with the correct experimentally measured intensity, while retaining the phase values and numerically propagate the distribution back to the image plane(4)We do not expect the calculated image intensity to be the same as the measured intensity distribution. Again we replace the numerically calculated distribution with the measured distribution(5)Steps (2)-(4) are repeated until the algorithm converges, i.e., when there is little difference between the measured and calculated intensity distributions. The similarity between the measured and numerically calculated distributions can be estimated with a root mean square error measure [46].

In [41], Fienup directly compares the behavior of this algorithm to gradient descent optimization approaches. The behavior and properties of these algorithms have been of significant interest in the field. It is generally said that for phase retrieval to work effectively one needs to have a good estimate of the initial phase. When that is the case the algorithms tend to converge to the nearest local minimum which in this case would be one of the global minimums. And hence the algorithm will return the real-physical phase distribution and not a solution that merely satisfies the numerical constraints of the problem but which has no physical meaning. In this manuscript we will explore this question of multiple phase solutions more closely, by examining a simplified version of the phase retrieval problem. A major outstanding issue with these iterative approaches is that algorithm often stagnates in a local minimum yielding a poor estimate for the phase. There is also the possibility that there are several different solutions to this ill-posed problem, at least in the 1D case [48].

In the iterative phase retrieval approach that we have just been discussing the intensity distributions are captured using CCD/CMOS arrays. And so in practice the intensity distribution is only recorded over a finite spatial extent and with a finite number of pixels which limits the spatial resolution and also introduces numerical sampling effects [37, 49]. The numerical procedures that are then used with these 2D arrays of intensity values tend to rely on fast Fourier transform (FFT) based algorithms for their implementation. This requires developing the software and numerical techniques for implementing these algorithms so that they are sufficiently fast and so that other sampling rules and numerical issues can be addressed [49]. This development, while necessary, tends to concentrate attention on the software implementation rather than on other physical and mathematical modeling aspects of the problem.

Therefore in this manuscript we deliberately take a different approach that emphasizes a theoretical analysis using a relatively simple model. Similar to the iterative phase retrieval approach we analyze how the phase may be estimated from two intensity measurements: one made in the image plane and one made in the Fourier plane. We do however make an important simplification. We assume that we can measure the intensity distribution in the Fourier plane over an infinite extent and ignore any sampling or discretization process that occurs due to the finite extent of the pixels in the CCD/CMOS array. Furthermore we assume that the intensity distribution in the image plane consists of a finite set of discrete point sources. With these two simplifications we can derive an analytical solution to the phase retrieval problem.

The manuscript is organized in the following manner: In Section 2, we show, under assumptions, which we have outlined in the text above and which are detailed in Section 2, how the phase values can in principle be determined using an analytical solution. To recover the phase values we must use an iterative procedure which has a finite but large number of steps. In Section 3, we examine an algorithm for implementing the iterative strategy derived in Section 2 and examine how the number of operations depends on the number of contributing point sources in the image plane. We show that there are multiple solutions to the phase retrieval problem. We then examine approaches for reducing the number of possible solutions. In Section 4, we finish with a brief discussion of the presented analysis and highlight some practical difficulties when extending the theoretical model to realistic real-world phase measurements.

#### 2. Analysis of an Idealized Experiment

We begin our analysis by examining the idealized optical system depicted in Figure 1. There is a set of distinct point sources, each with a known intensity value and hence a known magnitude, , where is the magnitude of point source . The phase value, , associated with each point source however remains the unknown physical parameter that we wish to measure. So to reiterate, the problem as depicted in Figure 1 is that we are given a set of point sources with known amplitudes and unknown phase values. We further assume that these points sources are uniformly spaced a distance from each other. We are also given a second measurement, the intensity distribution in the Fourier plane. What we have just described is very similar to the practical iterative phase retrieval experiment described in the Introduction. In a real experiment however these intensity distributions are recorded by CCD/CMOS arrays which have a finite extent, a finite number of pixels, and which sample of the intensity distribution they measure. These are critical factors when considering the performance limitations of these systems; see [33, 36]. They also complicate the theoretical analysis and so now we neglect these factors so that we can develop a simpler theoretical description of the phase retrieval problem.