## Abstract

The transport of intensity equation (TIE) is a two-dimensional second order elliptic partial differential equation that must be solved under appropriate boundary conditions. However, the boundary conditions are difficult to obtain in practice. The fast Fourier transform (FFT) based TIE solutions are widely adopted for its speed and simplicity. However, it implies periodic boundary conditions, which lead to significant boundary artifacts when the imposed assumption is violated. In this work, TIE phase retrieval is considered as an inhomogeneous Neumann boundary value problem with the boundary values experimentally measurable around a hard-edged aperture, without any assumption or prior knowledge about the test object and the setup. The analytic integral solution via Green’s function is given, as well as a fast numerical implementation for a rectangular region using the discrete cosine transform. This approach is applicable for the case of non-uniform intensity distribution with no extra effort to extract the boundary values from the intensity derivative signals. Its efficiency and robustness have been verified by several numerical simulations even when the objects are complex and the intensity measurements are noisy. This method promises to be an effective fast TIE solver for quantitative phase imaging applications.

© 2014 Optical Society of America

## 1. Introduction

Phase is an important component of an optical wavefield providing information of the refractive index, optical thickness, and the topology of the specimen. Phase retrieval is a central problem in many areas of physics and optics since the phase of a wavefield is not accessible directly. In contrast, it is relatively easy to measure the intensity of the wavefield; therefore it is tempting to seek methods for phase reconstruction from intensity measurements. The commonly known transport of intensity equation (TIE) proposed by Teague [1] is a non-interferometric method for the optical phase retrieval. This equation provides a relationship between intensity and phase of a light wave in the near Fresnel regime. With two or more intensity images taken at distinct planes orthogonal to the optical axis, solving the TIE directly gives the phase distribution, in a non-iterative manner. In the past few decades, TIE has found a variety of applications in adaptive optics [2], X-ray diffraction [3], electron- microscopy [4] and optical quantitative phase imaging [5, 6].

TIE is a second order elliptic partial differential equation for the phase function, and solving this equation does not appear to be difficult. Mathematically, the TIE phase retrieval is essentially a boundary value problem – seeking the solution of the partial differential equation subject to certain boundary conditions. The well-posed nature of this problem is well acknowledged in literature - with proper boundary conditions, one may find a unique (apart from an arbitrary additive constant) solution to the TIE if the intensity distribution is strictly positive [7]. Several method have been proposed to solve the TIE for the unknown phase: the Green’s function method [1, 8], the multi-grid method [9, 10], the Zernike polynomial expansion method [11, 12], and the fast Fourier transform (FFT) method [9, 13–15]. However, one important practical problem pervading these TIE solvers is that the associated boundary conditions are difficult to measure or to know a priori. For example, in Teague’s Green’s function solution [1], one needs to know the phase function at the region boundary as the Dirichlet boundary conditions, which is rather difficult without resorting to other means. The use of Neumann boundary conditions, studied by Roddier [2, 16] and later by Woods and Greenaway [8] has the advantage that the boundary values can be found from the intensity measurements at the pupil boundary. However, the results make one major assumption - that the intensity distribution should be uniform within the region [2, 16]. Moreover, distinguishing the boundary signal from the interior intensity derivative has been known to cause serious difficulties [11, 17]. Gureyev and Nugent [12] proposed to solve the TIE under non-classical boundary conditions – assuming the intensity distribution is non-zero inside the domain but drops to zero at the boundary. They proved that under such condition, TIE has a unique solution up to an additive constant. They then, employed orthogonal polynomial expansion (Zernike polynomial expansion for the circular region [11, 12] and fast Fourier transform (FFT) for the rectangular region [12, 13]) to solve the TIE without imposing any separate boundary conditions. However, the soft-edged intensity profile is hard to realize experimentally, and one may face the difficulty on how to define the correct domain as the intensity tends to zero as it approaches the boundary.

Another important issue among these TIE solvers is that only very few of them are sufficiently fast for reasonable sized image. The FFT-based methods [9, 13–15] are probably the most popular TIE solvers especially in the field of quantitative phase imaging [18–23] for its simplicity, efficiency, and applicability to rectangular domain. Though Gureyev and Nugent [12, 13] claimed that no boundary conditions are required for getting a unique (up to an arbitrary additive constant) solution of the TIE, the FFT-based methods have been shown to be quite sensitive to the boundary artifacts, in particular when the actual phase distribution covers the boundary of the image [15, 24, 25]. The root of the problem lies in the fact that the FFT solver implying periodic boundary conditions due to the cyclic nature of the discrete Fourier transform, which is often in contradiction with the actual situations. Throughout literature [18–23], a common way to avoid such boundary error is to consider the object to be isolated and placed in the center of the image [24]. However, this configuration does not reflect the general experimental conditions, and is not practical when the measured object is larger than the camera field of view. An alternative FFT-based method proposed by Volkov *et al.* [15] extends the measured intensity by a mirror padding scheme, which provides homogenous Neumann boundary conditions rather than the periodic ones. However, it still does not fundamentally solve the boundary artifacts problem since the boundary conditions are formed purely by a mathematical trick which is generally not physically grounded.

In this paper, we will solve the TIE with experimentally measurable inhomogeneous Neumann boundary conditions, through an efficient numerical algorithm using the fast discrete cosine transform (DCT). This approach not only addresses the above mentioned difficulties such as obtaining boundary conditions under non-uniform intensity distribution, distinguishing the boundary signal from the interior intensity derivative, and defining the correct domain for the solution, but also retains the major advantages of the FFT-based methods – simple and fast for the rectangular domain. Rather than conventional differential equation methods where the TIE is discretized directly, we derive the analytic integral solution via the Green’s function method, and then the DCT is applied to calculate the integral solution numerically. All the efforts are made to simplify the numerical calculation process as well as the enforcement of the experimentally obtained Neumann boundary conditions. Its good performance is demonstrated by several numerical examples.

## 2. Transport of intensity equation, boundary conditions, and uniqueness of the solution

Assume that a paraxial beam is propagating along the $z$ axis, and let the complex amplitude of the object be $\sqrt{I(r)}\mathrm{exp}[ik\varphi (r)]$. The derivative of intensity along the beam propagation direction, $z$ contains phase information that can be retrieved via TIE [1]:

#### 2.1 Boundary conditions and uniqueness of the solution

The TIE is an elliptic partial differential equation for the phase function, $\varphi $. We assume the region governed by the TIE to be a general open and bounded domain $\Omega \subset {\mathbb{R}}^{2}$ with a piecewise smooth boundary $\partial \Omega $. The intensity distribution $I$ is a continuous, non-negative function defined on the enclosure $\overline{\Omega}$ (consisting of the domain $\Omega $ plus its boundary $\partial \Omega $) and is smooth and strictly positive in $\Omega $. The longitudinal intensity derivative $\partial I/\partial z$ is assumed to be a smooth function in $\Omega $. The phase $\varphi $ is expected to be single-valued and smooth in $\overline{\Omega}$. Before deriving a solution to the TIE, the issues concerning its solvability (well-posedness) and the uniqueness of the solution must be addressed. The classic theory of elliptic partial differential equation states that the uniqueness of the TIE solution subject to certain boundary conditions [26]. We shall consider two classes of possible boundary conditions that could be applied to the solution of TIE: (1) *Dirichlet boundary conditions*. The values of phase function $\varphi $ are specified on the domain boundary

*Neumann boundary conditions*. The product of $I$ and the normal derivative of $\varphi $, are specified on the domain boundaryHere $g$ is a smooth function on the boundary $\partial \Omega $, $\partial \varphi /\partial n$ is the outward normal derivative. For the case of Dirichlet boundary conditions, the solution to the TIE always exists and is unique [26]. The case of the Neumann boundary conditions demands special attention because a solution may or may not exist (depending on whether the compatibility condition (Eq. (5)) is satisfied [7, 27]). Furthermore, it is clear that $\varphi $ = constant satisfies the homogeneous problem with $-k\partial I/\partial z=0$ and $g=0$, therefore if $\varphi $ is a solution, then $\varphi +C$ is also a solution for any constant $C$. In other words, the solution of the TIE subject to Neumann boundary conditions (Eq. (4)) should be unique up to an arbitrary additive constant, assuming that the solution exists in the first place. This constant is not essential for the phase retrieval problem. The proof of the uniqueness theorem is straightforward: we can first assume that two qualified solutions have been exhibited and then prove they can only differ from one additive constant using the maximum principle for elliptic equations [28].

#### 2.2 Compatibility condition and energy conservation law

For the Neumann boundary problem, there is an important compatibility condition that must be satisfied in order that a solution exists. Starting from the TIE (Eq. (1)), we integrate both sides of the equation over $\Omega $ and apply the divergence theorem to get the compatibility condition

It is more interesting to examine the physical picture described by Eq. (5): actually, it can be thought as an expression of *local energy conservation law* - the loss of energy (intensity) inside the region arising from energy flow across the region boundary. If we extend $\Omega $ to the unbounded space, the contour integral vanishes and Eq. (5) becomes

*law of energy conservation for unlimited free space*. In essence, the energy conservation expressed by Eqs. (5) and (6) is a universal law of physics, and their validity relies only on the paraxial approximation.

#### 2.3 Solving the TIE without any boundary values

Since solving the TIE requires boundary conditions either the Dirichlet boundary value or the normal derivative of the phase (Neumann boundary conditions) to ensure the uniqueness of the solution, and both kinds of boundary data are not easy to acquire experimentally, researchers are looking for ways to solve the TIE without taking any boundary measurements. Coincidently, all the efforts aim to find some way to nullify the LHS of Eq. (5), making boundary conditions unnecessary

Generally, Eq. (7) cannot be really valid on a bounded region unless the contour integral (LHS of Eq. (5)) vanishes. Physically, it is equivalent to the special case that there is no overall energy transfer across the region boundary.One simple and common way to satisfy this is to let the test object be centrally isolated in the field of view. Then the zero phase changes at the boundary (corresponding to homogeneous Neumann boundary conditions ${I\partial \varphi /\partial n|}_{\partial \Omega}=0$) can be assumed thus the LHS of Eq. (5) vanishes. In fact, in this case one can define not only the Neumann boundary conditions, but any uniform Dirichlet boundary conditions (constant phase at the boundary ${\varphi |}_{\partial \Omega}=C$) or periodic ones (the phase at the boundary repeats cyclically) as well. Gureyev and Nugent [12] suggested another way to eliminate the need of boundary conditions (LHS of Eq. (5)) by considering the case that the intensity $I>0$ inside the domain $\Omega $ but strictly vanishes at the boundary (actually it corresponds to a special case of the homogeneous Neumann boundary conditions ${I\partial \varphi /\partial n|}_{\partial \Omega}=0$ though not mentioned in [12]). Alternatively, without any additional requirement about the test object and experimental conditions, Volkov *et al.* [15] proposed a pure mathematical trick to nullify the energy flow across the boundary through appropriate symmetrization of input images. However, it assumes there is no energy dissipation through the image boundary for any class of images and objects, which is generally not physically grounded. To summarize former discussion, *without boundary value measurements does not mean that we can solve the TIE without imposing any boundary conditions, or more exactly, we have to confine our test object or experimental configuration to some implicit boundary conditions*. Obviously, the assumption of the test object to be isolated is apparently limiting and does not reflect general cases. Furthermore, the assumption is impractical when the object is larger than the image field of view. On the other hand, intensity measurement satisfying Gureyev and Nugent’s conditions tends to be difficult as well. One needs to use a special-designed soft-edged aperture or use specific illumination with intensity vanishing at the boundary, $\partial \Omega $. Furthermore, defining the boundary of the domain tends to be very difficult in practice since *very-near* zero and *at* zero may seem almost the same but can be entirely different *vis-a-vis* the well-posedness of the TIE.

## 3. Solution to transport of intensity under non-uniform intensity

As suggested in last section, because of its physical origin, solving the TIE essentially boils down to boundary value problems – to find a solution in $\Omega $ with conditions imposed on its boundary, $\partial \Omega $. However, the major obstacle for solving the TIE is that the associated boundary conditions are difficult to measure or know a priori. Therefore, in this section, a simple approach to get the boundary conditions and construct an analytic solution to the TIE using Green’s function method is suggested.

#### 3.1 Boundary values generation with a hard-edged aperture

Let us consider a particular intensity distribution in the object plane ($z=0$): the intensity $I$ is smooth (but can be non-uniform) in $\Omega $ but suddenly vanishes with a step-like discontinuity at the boundary, $\partial \Omega $. Physically, it is equivalent to capturing intensity images through a hard-edged aperture at the object plane:

Let us put aside these complications for a while and assumes we already have the separated longitudinal intensity derivative inside $\Omega $ and the amplitude of the ${\delta}_{\partial \Omega}$ along the boundary line through some peculiar way. When applying these measured quantities to solve the TIE (us the interior signal as the source term of the TIE (LHS of Eq. (1)) and take the amplitude of delta-function term in Eq. (11) as the function $g$ in Eq. (4)), we must make sure the Neumann boundary value problem is well-posed and yields a unique solution to the problem. Integrating both sides of Eq. (11) over all space and using Eq. (10), we obtain

#### 3.2 Analytic solution to the TIE with the Green’s function approach

In this subsection, we derive a special analytic solution to the TIE with the Neumann boundary conditions found from the intensity derivative boundary based on Green’s function approach. The reason why we do not use numerical methods directly is that these methods generally require the intensity derivative inside the aperture and the Neumann boundary value (delta-function term at the aperture edge) being treated separately. More specifically, in order to get the Neumann boundary value, one has to find some appropriate ways to extract the amplitude of the delta-function around the aperture edge. However, as we mentioned previously, the Eq. (11) is an only an idealized representation which is not suitable for direct application. The delta function ${\delta}_{\partial \Omega}$ is only defined on infinitely small support and takes infinite values, reflecting the fact that a physical quantity cannot be measured along one line. Although in reality the boundary signal assumes finite values and can be measured by a camera, it occupies a finite area, bringing about a series of non-trivial issues like how to define the boundary region, and how to separate the boundary values from the interior signals [11, 17]. In the following, we will introduce a way to bypass these difficulties.

Firstly, as suggested by Teague [1], an auxiliary function $\psi $ was introduced, which satisfies

This physical meaning of the vector field $I\nabla \varphi $ can be interpreted as the Poynting vector [14, 29], which is characterized by the scalar potential $\psi $. With this auxiliary function $\psi $, the elliptic TIE can be converted down to a standard Poisson equation [9, 13, 23]with the associated Neumann boundary conditions unchanged*the experimentally measured boundary value information has already been encapsulated in the area integral*. It can be seen that once the Green’s function $G\left(r,{r}^{\prime}\right)$ is determined, the solution of the Poisson equation (Eq. (15)) can be obtained from Eq. (20) with single area integral via Green's function; and we emphasize the Green's function here is solely determined by the shape of the region $\Omega $, and independent on experimental data.

Once $\psi $ is obtained by the area integral (Eq. (20)), the phase function $\varphi $ can be reconstructed from its determined gradient $\nabla \varphi ={I}^{-1}\nabla \psi $ uniquely up to an additive constant. Mathematically, phase integration is a well-defined problem and can be simply achieved by a line-by-line accumulation procedure. However, the phase integration is ill-posed (and not well defined) as the vector field ${I}^{-1}\nabla \psi $ is rarely integrable (see Sec. 6.2 for details). A common approach is to find the scalar function whose gradients are closest to the determined vector field ${I}^{-1}\nabla \psi $ in a least squares sense. It can be shown by calculus of variations [31, 32], the solution to this least squares error problem is equivalent to that of finding the solution of the following Poisson equation

with homogeneous Neumann boundary conditionsThe Green's function solution of the Neumann boundary problem (Eqs. (21) and (22)) is## 4. Fast numerical implementation with use of discrete cosine transform

Evidently, the analytic solution to TIE (Eqs. (20) and (23)) is not very convenient for numerical implementation since the Green's function is a four dimensional function - the direct calculation is memory and computationally demanding. For an $M\times N$ pixel image, $G\left(r,{r}^{\prime}\right)$ will contain ${M}^{2}{N}^{2}$ elements, which may require thousands of gigabyte memory to store. Besides, the direct numerical implementation of the integral requires expensive operations of $O\left({M}^{2}{N}^{2}\right)$complexity. Therefore, in this section, we focus on developing an efficient numerical implementation of the solution. More specifically, we consider the case of a rectangular domain since most detectors (like the CCD camera) are of rectangular geometry, thus in order to maximally fill the camera sensor area, using a rectangular aperture is a natural choice. It is found that the Green's function under Neumann boundary conditions with a rectangular region can be represented as an infinite-series expansion in terms of Fourier cosine harmonics, leading to a fast and efficient numerical implementation of the TIE solution with use of fast DCT. It should be noted that in practice this approach can be easily extended to the domains with arbitrary shapes (e.g. the circular or annular regions) as long as the corresponding eigen-functions of Laplacian can be obtained.

#### 4.1 Eigen-functions expansion of the Green's function

Although the analytic representation of the Green's function of Poisson equation under Neumann boundary conditions for a rectangular domain $[0,a]\times [0,b]$ can be found in many mathematical textbooks [30, 33], here we prefer to expand it into an infinite series to gain some insight to the efficient numerical implementation. Consider the Green's function in Eq. (20) associated with the Neumann condition problem, which obeys Eq. (18), we can expand it in terms of eigen-functions of Laplacian [33]

#### 4.2 Numerical implementation with fast discrete cosine transform

In practice, the intensity distributions are always recorded by a pixelated detector - the continuous function $\partial I/\partial z$ is discretized into a finite $M\times N$ rectangular grid cells of equal size. For convenience, in this subsection we still use the $r=\left(x,y\right)$ to represent the image traverse coordinate, but keep in mind they can only assume integral values, i.e., $x=0,1,\mathrm{...},M-1;$$y=0,1,\mathrm{...},N-1$. The pixel coordinate $\left(x,y\right)$ corresponds the grid cell centered at the location $\left(\Delta x\left(0.5+x\right),\Delta y\left(0.5+y\right)\right)$ in the continuous real space within the rectangular domain $[0,a]\times [0,b]$, where $\Delta x=a/M$, $\Delta y=b/N$ is the cell interval (pixel size) in the $x$ and $y$ direction, respectively. Besides the intensity derivative $\partial I/\partial z$, the cosine harmonics basis functions need to be discretized as well. As shown in Appendix B, the set of functions $\mathrm{cos}\left[m\pi \left(x+0.5\right)/M\right]\mathrm{cos}\left[n\pi \left(y+0.5\right)/N\right]$ form an orthogonal basis over the $M\times N$ discrete space while satisfies the homogeneous Neumann condition on the boundary. Correspondingly, the discretized form of Eq. (28) can be represented as

So far, the only remaining problem is to calculate the source term $\nabla \cdot \left({I}^{-1}\nabla \psi \right)$ in Eq. (23). As $\psi $ can be expanded into DCT basis vectors, its derivative can be efficiently calculated through DCT-based differential operators (all the relevant operators can be found in Appendix C). Once the gradient field $\nabla \psi =\left({\partial}_{x-DCT}\psi /,{\partial}_{y-DCT}\psi \right)$ is obtained from Eqs. (53) and (55), the vector field ${I}^{-1}\nabla \psi =\left({I}^{-1}{\partial}_{x}\psi ,{I}^{-1}{\partial}_{y}\psi \right)$can be easily calculated by point-wise multiplication. Finally, the divergence $\nabla \cdot \left({I}^{-1}\nabla \psi \right)$ can be obtained through

Combining the above analysis, the complete numerical implementation of the TIE solution can be summarized as the following steps:

- Calculate the intensity derivative term$-k\partial I/\partial z$ from the measured intensity images.
- Transform$-k\partial I/\partial z$ into cosine frequency domain by DCT (Eq. (43)).
- Reweight the DCT coefficients by ${\lambda}_{m,n}^{-1}$(Eq. (27)).
- Transform the modified DCT coefficient back to spatial domain to get $\psi $ (Eq. (44)).
- Calculate ${I}^{-1}\nabla \psi $.
- Transform $\nabla \cdot \left({I}^{-1}\nabla \psi \right)$ into cosine frequency domain by DCT.
- Reweight the DCT coefficients by ${\lambda}_{m,n}^{-1}$.
- The phase can be determined up to an additive constant by transforming the modified DCT coefficient back to spatial domain. Steps (8)-(10) are analogous to steps (2)-(4).

It should be noted that all the source intensity images should be strictly defined on the rectangular-shaped area $\overline{\Omega}$, which includes both the aperture boundary and the region inside it. However, it is evident that the above steps have redundant operations: the auxiliary function $\psi $ is the intermediate step in the solution and therefore should not be express explicitly. Taking out all the redundancies from the above steps and adopting the simple notation we used in Appendix C, the final algorithm is given in Algorithm 1.

Further, when the intensity is uniform in the defined region, Algorithm 1 can be further simplified as:

In this case, the TIE can be directly simplified as a Poisson equation without the need of the auxiliary function $\psi $, which corresponds to the case discussed in [8,16]. However, our method involves only one forward and one inverse DCT, which substantially reduces the computation and memory demand compared with the classical Green's function method [8]. As shown in Appendices B and D, the DCT as well as all DCT-based derivative operators involved can be easily calculated through FFT with even extension scheme. These methods are quite helpful to simplify the numerical implementation of the proposed algorithms since FFT is available in almost all computer packages, but fast algorithms for DCT and some other transforms like Eq. (54) are not. The two alternatives to Algorithm 1 and 2 based on FFT and even extension are given in Algorithm 3 and 4, respectively.

The meaning of the bar (‘-’) and angle brackets notations (‘< >’) used here can be found in Appendix D. At the first glance, the Algorithms 3 and 4 seem to be very similar with the previous work using FFT and even symmetrization [15]. However, one should note the key differences: our algorithm only applies to the experimental data captured with the presence of an aperture, which means the axial derivative of intensity generally contains a delta-function term at the boundary. It was derived from the Green's function method with the inhomogeneous Neumann boundary value measured at the aperture edge without imposing any implicit or prior assumptions on test objects. We found the solution to the Poisson equation can be expressed by cosine series expansion and then we use DCT as a fast numerical method for the Green's function solution. The FFT and even symmetrization schemes merely serve as alternative ways to perform DCT and its related computations. While previous work [15] used even symmetrization to nullify the contour integral (corresponding to the second terms of RHS of Eq. (17)) by making a silent assumption for the test object to be even symmetry, which is usually unrealistic especially when the object extends out of the image boundary, leading to non-optimal or even erroneous results.

Before proceeding further, several points on practical implementation of our algorithm need to be clarified: First, ${\lambda}_{m,n}^{-1}$ is not defined at the zero point of the DCT space, therefore, we multiply by zero instead. Second, since the intensity $I$ appears in the denominator, one should avoid dividing too small values to guarantee the stability of solution by setting a certain threshold value (e.g. 1% of the maximum value of $I$). Finally, one should correctly define the domain for each captured intensity image before carrying out the algorithms. It is trivial when the small defocus distance requirement is fulfilled (Eq. (11) is valid) - one may simply choose the pixels within aperture region $\Omega $ as well as those around the aperture boundary $\partial \Omega $. However, when the defocus distance is not sufficiently small (increasing the defocus distance is a common way to improve the SNR of the intensity derivative estimates [35]), some boundary signal will either go into the aperture zone or escape outside from the aperture edge due to the diffraction effect (Eq. (11) will no longer well hold). In this case, defining the proper domain seems problematic, and we will examine this point in Section 5.4.

## 5. Simulations and comparisons

Simulations were performed to test the proposed method. We simulated three challenging situations for TIE phase retrieval. In these simulations, the complex field was defined on a 256 × 256 grid with 2μm × 2μm pixel size. The wavelength is 632.8nm and the object is assumed to be larger than the camera field of view (200 × 200 pixels) which represents a more general situation than the common assumption in the literatures that the object is isolated and placed in the center of the camera field of view. Angular spectrum numerical propagation was used to generate two defocused images ${I}_{+}$ and ${I}_{-}$ with the distance of $\Delta z=\pm 10$μm. The intensity derivative $\partial I/\partial z$is estimated by central finite difference

#### 5.1 Pure phase object with uniform intensity distribution

First we focus on the case that the test object is a pure phase function with unitary intensity. As shown in Fig. 1(a), the phase function is predefined as $\varphi (x,y)=10{x}^{2}-10{y}^{2}-0.7x+2y+0.\text{82}$ where the lateral dimensions $x$ and $y$ are limited within a range from −0.2 mm to 0.2 mm in the camera field of view. Since the object is a pure phase one, $I$ is constant when no aperture exists, and can be pulled out of the divergence operator [21, 36, 37] from the TIE (Eq. (1)), resulting in a standard Poisson equation

This Poisson equation states that the intensity derivative is proportional to the Laplacian of the phase. Since the simulated phase function has zero Laplacian, the difference between two defocus images ${I}_{+}$ and ${I}_{-}$ provided a*null signal*as expected (Fig. 1(b)), therefore reconstructing the phase using these intensity measurements can never be achieved in this situation [38, 39]. This once again reveals the dependence of the TIE solution on the boundary conditions. By introducing an aperture (aperture size 180 × 180 pixels) located at the center of the object plane, intensity derivative signal can (only) be observed around the aperture boundary (Fig. 1(c)), which is consistent with our theoretical analysis (Eq. (11)). With use of the proposed method (Algorithm 2), the phase function was retrieved accurately with the RMSE 0.91% (Figs. 1(d) and 1(e)). Since only one DCT and inverse DCT is needed to retrieve the phase, the algorithm can be performed extremely fast. The computation time in this example was only 16 ms using a 2.5 GHz laptop in the MATLAB environment. In contrast, the direct computation of the convolution integral (Eq. (20)) requires a large amount of memory to store the Green’s function, therefore limits the maximum size of the phase matrix allowed to only 50 × 50 pixels. Even for such low resolution data, the direct computation took almost 12 seconds.

#### 5.2 Phase retrieval with non-uniform intensity distribution

To test the proposed algorithm in the case of non-uniform intensity distribution, we further consider the same phase map but with the non-uniform intensity distribution $I(x,y)=\mathrm{exp}\left[-\left({x}^{2}+{y}^{2}\right)/2{\sigma}^{2}\right]$ with $\sigma $ = 0.18 mm (Fig. 2(a)). In this case, the intensity derivative signal obtained is non-zero in the absence of the aperture, but still quite weak, as shown in Fig. 2(b). When the aperture was introduced, boundary signal can be observed at the aperture edge, which is much stronger compared with the derivative signal inside the aperture (Fig. 2(c)). Again, the phase distribution was faithfully retrieved by the proposed method (Algorithm 1) (Fig. 2(d)), with the residual error plotted in Fig. 2(e). The RMSE of the reconstructed phase is 2.23%, which is slightly higher than the uniform intensity case. The increased error may result from the “Teague’s assumption” (see Section 6.2 for details). The runtime of the algorithm in this example was 66 ms, which is nearly four times as the uniform intensity case (Algorithm 2). However, if one simply assumes constant intensity and applies Algorithm 2 to the same dataset when the intensity is in fact non-uniform, the resultant phase error will be much higher, as shown in Fig. 2(f). This further verifies that the intensity variations cannot be neglected even its distribution is relatively smooth.

To compare the proposed method with previously developed techniques, we applied the FFT-based
method developed by Paganin and Nugent [14] and the
FFT-based symmetrization method developed by Volkov *et.al* [15] to the same dataset but without the presence of the
aperture. Besides, for the FFT-based method, we also tested it with the same dataset as our
DCT-based algorithm in the presence of the aperture. The RMSEs of different methods are
summarized in the first row (Dataset 1) of Table
1.It can be seen that *only* the proposed algorithm reconstructed the
phase reliably, while all other methods failed to retrieve the phase function with zero
Laplacian correctly, leading to very large RMSE values.

#### 5.3 Complex phase and intensity distribution

The next simulation focuses on more complex phase map and intensity distribution. We
consider the phase profile of letters ‘TIE’ (Fig.
3(a)) with the non-uniform intensity distribution of letters ‘COLE’
(Fig. 3(b)). Figures
3(c) and 3(d) plot the longitudinal intensity
derivative distribution without/with the presence of the aperture, respectively. The two
derivative images look quite similar at first glance. However, a closer inspection on the
enlarged area (Figs. 3(e) and 3(f)) reveals that *the aperture generated additional boundary
signals at the aperture edge which cannot be observed for the no aperture case*.

Figure 4 gives the pictorial comparison of different solutions and the RMSE values are presented in the second row (Dataset 2) of Table 1. It can be seen that the conventional FFT method could not retrieve the phase correctly no matter with (Figs. 4(b) and 4(f)) or without (Figs. 4(c) and 4(d)) the aperture since it simply assumes periodic boundary conditions, which is of course not true in this situation. Furthermore, the boundary error of the even symmetrization method was also significant in the absence of the aperture (Figs. 4(d) and 4(h)). In contrast, the proposed method provided a faithful result with much higher accuracy (RMSE 2.75%), proving that boundary errors can only be eliminated by the proposed DCT-based method with the practical boundary value information obtained on the aperture edges.

#### 5.4 Phase retrieval under noise and large defocus distance

In the last simulation, we tested the accuracy and stability of the proposed method in the presence of noise. Different levels of pseudo-random Gaussian noise with standard deviations (STD) of 0.0001, 0.001 and 0.01 were generated and superimposed on each intensity image of the data set from Section 5.3. To clearly demonstrate the noise levels, we show the calculated noisy intensity derivatives in Figs. 5(a)–5(c) (the defocus distance is $\Delta z=\pm 10$μm, same as the previous simulations). The reconstructed phases and corresponding residual errors are shown in Figs. 5(d)–5(f) and Figs. 5(g)–5(i), respectively. The RMSE values of the reconstructed phases are further given in the first row (defocus distances 10μm) of Table 2It can be seen that when the noise level was low, the current method could still recover the phase fairly accurately, with a slight increase in RMSE, from 2.23% to 3.14%. However, with increase of noise level, low-frequency artifacts can be clearly observed from the recovered phases. Especially when the noise STD reached 0.01, the original phase signal was almost submerged by the cloud-like artifacts, leading to significant error (RMSE 213.44%).

The high sensitivity to low-frequency phase error is a common and notorious problem in the TIE because the inverse Laplacian operator is ill-conditioned at near zero-frequency [35]. The common countermeasures include regularization treatments [4], using more intensity measurements [21, 22], or simply increasing the defocus distance [35]. By increasing the defocus distance $\Delta z$ from 10μm to 100μm, the derivative SNR is almost increased by 10 times since the image noise is amplified by the inverse of defocus distance $\Delta z$ in the derivative estimates, as verified by comparing the two intensity derivative distributions shown in Figs. 5(b) and 6(a).However, on the other hand, increasing $\Delta z$ also compromised the linearity assumption in the finite-difference derivative approximation (Eq. (32)), leading to non-linearity errors [21, 22]. From the enlarged region (Fig. 6(b)), we can clearly see the estimated derivative became blurred compared with the more accurate one obtained with $\Delta z$ = 10μm (Fig. 3(f)). Moreover, some boundary signals escaped from the aperture edge and went into the exterior region, suggesting Eq. (11) is no longer well hold. In this case, a new integration region ${\Omega}^{+}$ should be defined to make sure all the escaped signals can be enclosed in to fulfill the energy conservation law (Eq. (13)). Obviously, the choice of region ${\Omega}^{+}$ satisfying Eq. (13) is not unique and one can even extend ${\Omega}^{+}$ to the unbounded space. However, since the Green's function $G\left(r,{r}^{\prime}\right)$ in the integral (Eqs. (20) and (23)) is only defined on the aperture region; we should choose ${\Omega}^{+}$ as close to $\overline{\Omega}$ as possible to minimize this region mismatch. Therefore, a simple searching process was used to quickly define the appropriate integration region ${\Omega}^{+}$ by checking the integral value (Eq. (13)) starting from the initial region ${\Omega}^{+}=\overline{\Omega}$ and then gradually expanded outwards until it is sufficiently close to zero. The new integral region ${\Omega}^{+}$ (190 × 190 pixels) found by this procedure is indicated by the solid black line in Fig. 6(b), which is slightly larger than the original aperture region (180 × 180 pixels). Since the intensity $I\left(r\right)$ is defined only on the closure $\overline{\Omega}$, the gap (intensity zeros) region between $\overline{\Omega}$ and ${\Omega}^{+}$can duplicate the intensity values of the outermost pixels within the aperture region (or simply restrict the domain of phase integration in Eq. (23) to $\overline{\Omega}$ only). The reconstructed phase and the residual error are shown in Fig. 6(c) and 6(d), respectively, which demonstrates a great improvement in phase accuracy compared with the result obtained when $\Delta z$ = 10μm (Fig. 5(e)), with the RMSE value decreasing from 29.38% to 9.30%. But for the high-noise-level condition (noise STD 0.01), we need to further increase the defocus distance $\Delta z$ to 500μm. In this case, although the SNR of the derivative signal was improved greatly compared with Fig. 5(c), the nonlinear effect was significant – the estimated derivative became quite blurry, as shown from the enlarged region (Fig. 6(f)). It is obvious that we cannot expect the result to be very accurate for such a situation. Followed by a similar procedure, we found the integration region ${\Omega}^{\prime}$ in this case should be extended to 198 × 198 pixels indicated by the solid black line in Fig. 6(f). Once again, the reconstructed phase demonstrates a significant reduction in low-frequency artifacts, with the RMSE value dropping from 213.44% to 17.92%. However, the blurring effects can be easily perceived in the reconstructed phase and residual error map, which means the non-linearity error was quite dominant and the $\partial I/\partial z$ obtained by Eq. (32) was no longer a valid estimate of axial intensity derivative. However, even under such conditions, the RMSE of our algorithm is still much lower than all the other algorithms under ideal, noise-free conditions. These results confirm the stability of our algorithm with respect to the noise, and its validity under large defocus distances and nonlinearity errors.

## 6. Discussions

#### 6.1 Relation between FFT-based method

The main idea of the present method is to consider the solution of the TIE as an inhomogeneous Neumann boundary value problem - with the Neumann boundary value experimentally measurable around a hard-edged aperture. The validity of boundary value is guaranteed by the energy conservation law so that our method represent a more general case and can be applied to any complex object without imposing any implicit or prior assumptions. Conversely, as we discussed in Sec. 2.3, if we are able to make some safe assumptions on the test object or experimental configuration such as constant phase at the boundary (which are a quite common assumptions to ensure the validity the FFT-based methods), the boundary measurements and the aperture become unnecessary because there would be no boundary signal at the aperture edge (see Eq. (9)). In this case, one can define not only the Neumann boundary, but any uniform Dirichlet boundary conditions or periodic ones as well. The eigen-functions of the Laplacian can then be more general Fourier harmonics [33] (can be cosine, sine, or exponential harmonics), and thus there will be no essential difference between the proposed DCT-based solution and the traditional FFT-based solution.

#### 6.2 Inaccuracy results from “Teague’s assumption”

It has to be mentioned that the TIE solution presented here is based on the “Teague’s assumption” that there exists an auxiliary function $\psi $ satisfying Eq. (14). However, as some previous work suggested, Teague’s assumption may not be always valid [9, 40]. According to the Helmholtz’s theorem, the vector field $I\nabla \varphi $can be decomposed as

where $rot\eta $ is the vector field $rot\eta =\left(\partial \eta /\partial y,\partial \eta /\partial x\right)$ and $\eta $ is a scalar function. Compared with Eq. (14), it can be seen that the rotational term $rot\eta $ was simply ignored in Teague’s assumption [9], postulating that $I\nabla \varphi $ must be potential and irrotational. A detailed discussion about the validity and the inaccuracy resulting from “Teague’s assumption” can be found in [40]. Using the identity $\nabla \cdot rot\eta =0$, the first Poisson equation (Eq. (15)) is still valid and this approximation does not affect accurate reconstruction of function $\psi $. However, the missing rotational term does affect the solution of the second Poisson equation (Eq. (21))## 7. Conclusions

In conclusion, we have presented a method for TIE phase retrieval with use of DCT. The new method has several substantial advantages: First, it clearly defines the TIE as an inhomogeneous Neumann boundary value problem with the boundary value experimentally measurable around a hard-edged aperture. Second, the method is valid for the case of non-uniform intensity distribution. Third, the measured intensity data are treated as one entity, without requiring special-purpose detection scheme to separate boundary values from the interior signals. Finally, for rectangular region, the DCT based numerical implementation is quite simple, fast, and of low memory usage. Its good performance has been confirmed by several numerical simulations even when the objects are located at the image borders. Those precise results are clearly essential when dealing with quantitative phase imaging, which was the main driver for this research. The method has already been tested experimentally with the application to micro-optics characterization and found to be robust to many source of experimental error. The experimental results will be published in a subsequent paper.

## Appendix A: Cosine series expansion and discrete cosine transform

The preceding Appendixes served to motivate the development of DCT as an efficient tool for solving the two Poisson equations, i.e., calculating the two area integrals of Eqs. (20) and (23) related to the TIE solution. Considering a band limited continuous signal $f(x)$ defined on $[0,a]$, it can be expended into an infinite but convergent series of cosine basis functions

The coefficient for this expansion can be represented as*N*-dimensional vector space and meanwhile fulfills the Neumann at the region boundary.

Similar to the cosine transform of one dimension, for the two-dimensional (2-D) data sequence$\left\{f(x,y),x=0,1,\mathrm{...},N-1,y=0,1,\mathrm{...},M-1\right\}$, their DCT and inverse transform are defined respectively as:

## Appendix B: Calculate DCT based on FFT

The direct implementation of DCT involves four nested loops with complexity $O({M}^{2}{N}^{2})$ for an $M\times N$ array, which is prohibitive in most applications. As a result, various fast and efficient DCT algorithms have been devised and studied extensively [43]. However, here we only discuss one fast algorithm that utilizes the FFT to compute the DCT and its inverse [44] due to its simple implementation and relevance to one previous work [15]. The fast computation procedure consists of extending the input sequence of *N* samples to an 2*N*-point sequence with even symmetry, taking an 2*N*-point DFT, and saving *N* terms from it. The even extension is defined as

*N*-point DFT vector of $\overline{f}(n)$ is

*N*-point DFT coefficients $\overline{F}(k)$ from the DCT coefficients $F(k)$ as follows:

*N*-point inverse DFT, the reconstructed data by inverse DCT can be obtained in the first

*N*elements of the inverse DFT result.

## Appendix C: Derivative calculation based on DCT

In this section, expressions are derived for computing derivative operations on a signal, from the DCT coefficients of the original signal’s samples. Applying the derivative operation to both sides of Eq. (36), and adopting the notation ${\Psi}_{k}^{(l)}\left(x\right)$ mean the *l*th order derivative of the basis functions of the transform. Then *l*th order derivative of $f$can be deduced as

## Appendix D: Calculate DCT-based derivatives using FFT

In appendix C, we present several formulations for differentiation differential operators based on DCT. It has been found that for operators related to only even order derivatives, such as ${\partial}_{x-DCT}^{2}f$ and ${\nabla}_{{}^{DCT}}^{{}^{-2}}f$, the reconstruction is a simple inverse DCT. However, this rule does not work for odd order derivatives. The direct calculation of Eqs. (55) and (56) requires expensive computation $O({M}^{2}{N}^{2})$ for an $M\times N$ matrix. In Appendix B, we demonstrate a handy fast algorithm to compute DCT based on 2*N* point FFT. Accordingly, here we present the corresponding fast algorithms to calculate DCT-based differential operators using FFT.

We again consider the 2*N*-point sequence $\overline{f}(n)$ extended from the extending the input *N* samples according to Eq. (45) with Fourier coefficients $\overline{F}(k),\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}k=0,1,\mathrm{...},2N-1$ defined in Eq. (46). Using the relationship between DCT coefficients and the 2*N*-point DFT coefficients (Eqs. (47) and (48)), we can deduce that if we modify the $\overline{F}(k)$ as the following rules,

*N*-point inverse FFT, then the first

*N*elements of the reconstruction result corresponds to the sampled derivative calculated by DCT-based method (Eq. (51)).

*N*elements in the reconstruction result by FFT corresponds to the sampled second order derivative calculated by DCT-based method (Eq. (52)).

## Acknowledgments

We are grateful to the anonymous reviewers for their constructive suggestions. This project was supported by the Research Fund for the Doctoral Program of Ministry of Education of China (No. 20123219110016), the National Natural Science Foundation of China (No. 61271332). C. Zuo acknowledges the financial support from China Scholarship Council (No. 201206840009).

## References and links

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

**2. **F. Roddier, “Curvature sensing and compensation: a new concept in adaptive optics,” Appl. Opt. **27**(7), 1223–1225 (1988). [CrossRef] [PubMed]

**3. **K. A. Nugent, T. E. Gureyev, D. J. Cookson, D. Paganin, and Z. Barnea, “Quantitative phase imaging using hard X rays,” Phys. Rev. Lett. **77**(14), 2961–2964 (1996). [CrossRef] [PubMed]

**4. **K. Ishizuka and B. Allman, “Phase measurement of atomic resolution image using transport of intensity equation,” J. Electron Microsc. **54**(3), 191–197 (2005). [CrossRef] [PubMed]

**5. **N. Streibl, “Phase imaging by the transport equation of intensity,” Opt. Commun. **49**(1), 6–10 (1984). [CrossRef]

**6. **C. Zuo, Q. Chen, W. Qu, and A. Asundi, “Noninterferometric single-shot quantitative phase microscopy,” Opt. Lett. **38**(18), 3538–3541 (2013). [CrossRef] [PubMed]

**7. **T. E. Gureyev, A. Roberts, and K. A. Nugent, “Partially coherent fields, the transport-of-intensity equation, and phase uniqueness,” J. Opt. Soc. Am. A **12**(9), 1942–1946 (1995). [CrossRef]

**8. **S. C. Woods and A. H. Greenaway, “Wave-front sensing by use of a Green’s function solution to the intensity transport equation,” J. Opt. Soc. Am. A **20**(3), 508–512 (2003). [CrossRef] [PubMed]

**9. **L. J. Allen and M. P. Oxley, “Phase retrieval from series of images obtained by defocus variation,” Opt. Commun. **199**(1–4), 65–75 (2001). [CrossRef]

**10. **S. V. Pinhasi, R. Alimi, L. Perelmutter, and S. Eliezer, “Topography retrieval using different solutions of the transport intensity equation,” J. Opt. Soc. Am. A **27**(10), 2285–2292 (2010). [CrossRef] [PubMed]

**11. **T. Gureyev, A. Roberts, and K. Nugent, “Phase retrieval with the transport-of-intensity equation: matrix solution with use of Zernike polynomials,” J. Opt. Soc. Am. A **12**(9), 1932–1942 (1995). [CrossRef]

**12. **T. E. Gureyev and K. A. Nugent, “Phase retrieval with the transport-of-intensity equation. II. Orthogonal series solution for nonuniform illumination,” J. Opt. Soc. Am. A **13**(8), 1670–1682 (1996). [CrossRef]

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

**14. **D. Paganin and K. A. Nugent, “Noninterferometric Phase Imaging with Partially Coherent Light,” Phys. Rev. Lett. **80**(12), 2586–2589 (1998). [CrossRef]

**15. **V. V. Volkov, Y. Zhu, and M. De Graef, “A new symmetrized solution for phase retrieval using the transport of intensity equation,” Micron **33**(5), 411–416 (2002). [CrossRef] [PubMed]

**16. **F. Roddier, “Wavefront sensing and the irradiance transport equation,” Appl. Opt. **29**(10), 1402–1403 (1990). [CrossRef] [PubMed]

**17. **I. W. Han, “New method for estimating wavefront from curvature signal by curve fitting,” Opt. Eng. **34**(4), 1232–1237 (1995). [CrossRef]

**18. **A. Barty, K. A. Nugent, D. Paganin, and A. Roberts, “Quantitative optical phase microscopy,” Opt. Lett. **23**(11), 817–819 (1998). [CrossRef] [PubMed]

**19. **E. D. Barone-Nugent, A. Barty, and K. A. Nugent, “Quantitative phase-amplitude microscopy I: optical microscopy,” J. Microsc. **206**(3), 194–203 (2002). [CrossRef] [PubMed]

**20. **S. S. Kou, L. Waller, G. Barbastathis, and C. J. R. Sheppard, “Transport-of-intensity approach to differential interference contrast (TI-DIC) microscopy for quantitative phase imaging,” Opt. Lett. **35**(3), 447–449 (2010). [CrossRef] [PubMed]

**21. **L. Waller, L. Tian, and G. Barbastathis, “Transport of intensity phase-amplitude imaging with higher order intensity derivatives,” Opt. Express **18**(12), 12552–12561 (2010). [CrossRef] [PubMed]

**22. **C. Zuo, Q. Chen, Y. Yu, and A. Asundi, “Transport-of-intensity phase imaging using Savitzky-Golay differentiation filter--theory and applications,” Opt. Express **21**(5), 5346–5362 (2013). [CrossRef] [PubMed]

**23. **C. Zuo, Q. Chen, W. Qu, and A. Asundi, “High-speed transport-of-intensity phase microscopy with an electrically tunable lens,” Opt. Express **21**(20), 24060–24075 (2013). [CrossRef] [PubMed]

**24. **J. Martinez-Carranza, K. Falaggis, T. Kozacki, and M. Kujawinska, “Effect of imposed boundary conditions on the accuracy of transport of intensity equation based solvers,” Proc. SPIE **8789**, 87890N (2013).

**25. **J. Frank, S. Altmeyer, and G. Wernicke, “Non-interferometric, non-iterative phase retrieval by Green’s functions,” J. Opt. Soc. Am. A **27**(10), 2244–2251 (2010). [CrossRef] [PubMed]

**26. **D. A. Gilbarg and N. S. Trudinger, *Elliptic Partial Differential Equations of Second Order* (Springer, 2001), Vol. 224.

**27. **K. Bhamra, *Partial Differential Equations: An Introductory Treatment With Applications* (PHI, 2010).

**28. **R. Courant and D. Hilbert, “Potential theory and elliptic differential equations,” in *Methods of Mathematical Physics* (Wiley-VCH Verlag GmbH, 2008), pp. 240–406.

**29. **C. Zuo, Q. Chen, and A. Asundi, “Light field moment imaging: comment,” Opt. Lett. **39**(3), 654 (2014). [CrossRef] [PubMed]

**30. **D. W. Trim, *Applied Partial Differential Equations* (PWS-Kent, 1990).

**31. **A. Agrawal, R. Raskar, and R. Chellappa, “What is the range of surface reconstructions from a gradient field?” in *Computer Vision–ECCV 2006* (Springer, 2006), pp. 578–591.

**32. **A. Talmi and E. N. Ribak, “Wavefront reconstruction from its gradients,” J. Opt. Soc. Am. A **23**(2), 288–297 (2006). [CrossRef] [PubMed]

**33. **F. Morse, *Methods of Theoretical Physics* (1981), Vols. 1 and 2.

**34. **W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, *Numerical Recipes in Fortran 77: The Art of Scientific Computing* (Cambridge University, 1992).

**35. **D. Paganin, A. Barty, P. J. McMahon, and K. A. Nugent, “Quantitative phase-amplitude microscopy. III. The effects of noise,” J. Microsc. **214**(1), 51–61 (2004). [CrossRef] [PubMed]

**36. **C. Dorrer and J. D. Zuegel, “Optical testing using the transport-of-intensity equation,” Opt. Express **15**(12), 7165–7175 (2007). [CrossRef] [PubMed]

**37. **L. Tian, J. C. Petruccelli, and G. Barbastathis, “Nonlinear diffusion regularization for transport of intensity phase imaging,” Opt. Lett. **37**(19), 4131–4133 (2012). [CrossRef] [PubMed]

**38. **C. Campbell, “Wave-front sensing by use of a Green’s function solution to the intensity transport equation: comment,” J. Opt. Soc. Am. A **24**(8), 2480–2482 (2007). [CrossRef] [PubMed]

**39. **S. C. Woods, H. I. Campbell, and A. H. Greenaway, “Wave-front sensing by use of a Green's function solution to the intensity transport equation: reply to comment,” J. Opt. Soc. Am. A **24**(8), 2482–2484 (2007). [CrossRef]

**40. **J. A. Schmalz, T. E. Gureyev, D. M. Paganin, and K. M. Pavlov, “Phase retrieval using radiation and matter-wave fields: Validity of Teague's method for solution of the transport-of-intensity equation,” Phys. Rev. A **84**(2), 023808 (2011). [CrossRef]

**41. **V. V. Volkov and Y. Zhu, “Lorentz phase microscopy of magnetic materials,” Ultramicroscopy **98**(2-4), 271–281 (2004). [CrossRef] [PubMed]

**42. **N. Ahmed, T. Natarajan, and K. R. Rao, “Discrete cosine transform,” IEEE Trans. Comput. **C-23**(1), 90–93 (1974). [CrossRef]

**43. **E. Feig and S. Winograd, “Fast algorithms for the discrete cosine transform,” IEEE Trans. Signal Process. **40**(9), 2174–2193 (1992). [CrossRef]

**44. **J. Tribolet and R. E. Crochiere, “Frequency domain coding of speech,” IEEE Trans. Acoust. Speech Signal Process. **27**(5), 512–530 (1979). [CrossRef]

**45. **R. Reeves and K. Kubik, “Shift, scaling and derivative properties for the discrete cosine transform,” Signal Process. **86**(7), 1597–1603 (2006). [CrossRef]

**46. **R. N. Bracewell, *The Fourier Transform and Its Applications* (McGraw-Hill, 1986), Vol. 3.