1. Introduction

In most optical systems, the measured information is usually the intensity of the image (signal) or its diffraction pattern. Hence, the substantial information encoded in the phase vanishes. Reconstructing both the magnitude and the phase information from the measured intensity in an optical field has been an important widely investigated issue [1-8]. Most of the previous methods were based on the two measured intensities or diffraction patterns in the Fourier transform domain. For example, Gerchberg and Saxton [2] proposed an iterative phase-retrieval algorithm that bounces back and forth between two planes that have the Fourier transform relation to each other. Cong et al. [9] proposed a different recursive algorithm to perform the phase retrieval from two intensities in the Fresnel transform system. This method can offer a multifarious choice of two intensities at arbitrary locations in the freespace diffraction planes and minimize the hardware requirement to implement it without the lenses needed in the Fourier domain. Also the oversampling phasing method was proposed and has been successfully applied to the coherent X-ray microscopy for 2D or 3D images of nanocrystal or noncrystalline specimens [10-12]. However, all aforementioned algorithms still need two intensities even in different transform systems. For this reason, a simpler and more efficient algorithm that needs only one single Fresnel transform intensity to reconstruct the signal (both the phase and the magnitude information) is presented. In the following analysis we consider, for simplification, the case of one dimensional function and one-dimensional transform. The generalization to two dimensional separable functions is straightforward.

2. Theory

The Fresnel transform (FrT^z) of a function f(x_o) is given by

where z is the propagation distance, λ is the wavelength, x_o is the variable of the coordinate of the input plane, and x_p is the variable of the coordinate of the output plane. As illustrated in Fig. 1, the optical field in the output plane x_p y_p is F^z (x_p), but only the intensity |F^z(x_p)|² is recorded. The aim is to rebuild the signal f(x_o) (both the phase and the magnitude information) from only one |F^z(x_p)|².

 

Fig. 1. The configuration to illustrate the Fresnel transform.

Download Full Size | PDF

In most of the realistic optical systems, the input and output signals can be regarded as being approximately both space limited and band limited. Assume the signal extension is 2Δx_o and there are N sample points with interval δx_o = 2Δx_o/N. This must satisfy the Nyquist sample theory that δx_o is less than the inverse of twice the signal’s bandwidth because the phase information is uniquely embedded inside the diffraction pattern when the pattern is sampled at a spacing finer than the Nquist frequency (i.e. oversampling) [13-15]. The sequence {f(mδx ₀) : m = N/2 ,…, N/2 -1 } can be generated from the continuous f(x _o). Consequently, we define the discrete Fresnel transform (DFrT^z) of f(mδx_o) as

where δx_o and δx_p are the sampling period in x_o and x_p space, m and n are integers, and κ(m, n, z) is the kernel of the DFrT^z as

The inverse DFrT^z (IDFrT^z) is given by

where δx_p = λz/(Nδx_o), and the superscript ^* denotes the complex conjugation.

The correlation property for DFrT^z can be represented, from Eqs. (2)-(4), as follows [9]

where f(m) and F^z(n) denote f(mδx_o) and f^z(nδx_p), respectively.

We now consider the problem of reconstructing a function f(m) from only the single intensity of its discrete Fresnel transform |F^z(n)|².

In a discrete transform system, the zero padding skill is employed to avoid the time aliasing which results from the periodic (circular) correlation. In order to obtain the correct linear correlation relation, we set the zero padding form of f(m) as

where the original data are still stored in k = -N/2 ,…, N/2 -1. With the help of Eqs. (5) and (6), the correlation property can be rewritten as

where δx_p = λz/(2Nδx_o), and let R(k) be the entire right-hand side of Eq. (7) as

If the autocorrelation lag is set as k = N - 1 , from Eqs. (7) and (8), we have

Next, for the lag k = N - 2,

Similarly, we derive the relation, for the lags k = N - m for m≥3, as

We can sort and arrange all the product factors without the phase term for different lags of k = N - m as in part A below. Each row denotes the same lag. For example, the first row is for k = N - 1; the second row is for k = N - 2 , and so on. The symbol A(i,j) (the factor in ith row, jth term) is used to help identifying the factor; for example, A(2,1) is the product factor f ^*(-N/2)f(N/2 - 2) because it is on the second row and the first term.

It is obvious that for all product factors, only f ^* (-N/2)f(N/2-1) [=A(1,1)] in the first row can be obtained directly from Eq. (9); thus the other factors (not to mention the sequence{f(m) : m = -N/2 , …, N/2-1}) can not be found using Eqs. (10) and (11). This is why the previous study claimed that the phase retrieval algorithm needs two intensities or diffraction patterns to reconstruct the signal [9].

To create more relations for the product terms to reconstruct the signal sequence, a property called the complex-convolution is developed as

The above equation can be obtained from Eq. (7) and using the symmetrical property of the signal (f(m) = f(-m)). Setting k = -N in Eq. (12), and let R′(k) be the entire right-hand side of Eq. (12) as

We have the following equation as

Next, for k = -N+1,

For the lags k = -N+m with m≥2 , the generalized expression is

As we do in part A, from Eqs.(14)-(16), the complex-convolution factors without phase term can be sorted and arranged for different lags k = -N + m as in part B below. Each row denotes the same lag, such as the first row for k = -N , second row for k = -N+1 , and so on; and the symbol B(i,j) is used to identify the indicated term. Again, only B(1,1) (=f ^*(-N/2)f(-N/2)) in the first row can be obtained directly by Eq.(14), but these product factors are very important for our aim as explained later.

3. The recursive algorithm

Now, we introduce the signal reconstruction algorithm and the procedures are explained as in the following. Step 1: The exact values of A(1,1) (= f ^* (-N/2)f(N/2-1)) and B(1,1) (= f ^* (-N/2)f(-N/2)) are calculated from Eqs. (9) and (14). It is found that B(2,1) (= f ^* (-N/2)f(-N/2+1)) is equal to A(1,1) (= f ^* (-N/2)f(N/2-1)) because f (m) is symmetrical and the condition f (-N/2 +1) = f (N/2-1) is satisfied. Step 2: It is also noted that B(2,2) = B(2,1)^*; therefore the complex conjugate operation can be used to gain the factor B(2,2) as for f ^* (-N/2 +1)f(-N/2) =[f ^*(-N/2)f(-N/2+1)]^*. After that, the factor B(3,2) (= f ^* (-N/2+1)f(-N/2+1)) in third row of part B can be calculated using the first two rows as (B(3,2) = B(2,1)×B(2,2)/B(1,1)), or explicitly as

With the help of the symmetrical condition f (-N/2+1) = f (N/2-1) , A(2,2) is equal to the B(3,2). Step 3: A(2,1) (= f ^* (-N/2)f(N/2-2)) can be obtained when A(2,2) (= f ^* (-N/2+1)f(N/2-1)) is substituted in Eq. (10). Once this is found, the factors B(3,1) (= f ^* (-N/2)f(-N/2+2)) and B(3,3) (= f ^* (-N/2+2)f(-N/2)) in part B can be obtained immediately for the symmetrical and conjugation properties (A(2,1)=B(3,1)= B(3,3)^*) as mentioned above. The factors in the third row of part B are known and then continue the same procedure to find the factors B(4,2) (= f ^* (-N/2+1)f(-N/2+2)) and B(4,3) (= f ^* (-N/2+2)f(-N/2+1)) in the forth row of part B as B(4,2) = B(3,1)×B(3,2)/B(2,1) and B(4,3) = B(3,2)×B(3,3)/B(2,2). Return these two values into the third row of part A (A(3,2) = B(4,2), A(3,3) = B(4,3)) and utilize Eq. (11) for k = N-3 to solve the factor A(3,1) (= f ^* (-N/2)f(N/2-3)) in the third row of part A because A(3,2) and A(3,3) are already known. This technique can be used repeatedly and the algorithm flow chart is depicted in Fig. 2. This recursive process continues until the value k = N-(N/2+1) is reached, which is the last row (the k = N/2-1 row as shown in part A triangle) and every product term in part A is found. Finally we can take the product of the two known terms f ^*(-N/2)f(0) and f ^* (0)f(N/2-1) (they are the first tem in k = N/2 row and the last term in k = N/2-1 row as shown in part A triangle respectively) and use Eq. (9) to get.

Therefore |f(0)| can be determined using the above equation, and only its phase is free. Choosing an arbitrary value for its phase does not have any observable effect on the recovered signal since it can be used as a reference. Now f(0) is known and refer to the first term and last term from bottom to head in part A. We can divide f ^*(-N/2)f(0) (the first term in k = N/2) by f(0) to obtain f ^*(-N/2) (thus f (-N/2)). And we can also divide f ^* (0)f(N/2-1) (the last term in k = N/2-1 row) by f ^* (0) to obtain f(N/2-1). For the same reason, the value of f(1) can be found when f ^* (-N/2)f(1) (the first term in k = N/2+1 row) is divided by f ^* (-N/2). Also the value of f ^* (-1) (thus f (-1)) is found when f ^* (1)f(N/2-1) (the last term in k = N/2 row) is divided by f (N/2-1)). This process is repeated and all the sequence [f (m) : m = -N/2, -N/2+2, ⋯, N/2-1] can be calculated. The signal is thus reconstructed.

 

Fig. 2. Signal reconstruction algorithm based on a single Fresnel transform intensity

Download Full Size | PDF

4. Numerical results

A symmetrical signal f(x_o) = exp[-16x ² _o+j 0.5sin(64πx ² _o)] is employed as an example for demonstrating the practicability of the proposed algorithm. The conditions for wavelength λ = 632.8 nm and the propagation distance z = 0.5 m are taken and the number of sampling points is N = 64. Only one intensity |F^z(x_p)|² is used in our algorithm to reconstruct f(x_o). The recovered magnitude and phase information are indicated in Figs. 3(a) and 3(b), respectively. In Fig. 3, the solid curves correspond to the original signal f(x_o). The open circles indicate signal recovered by our recursive method. From the results, it is evident that this algorithm performs very well and the reconstructed signal is almost the same as the original. This algorithm is also applied to another more complicated two-dimensional signal with the form f(x_o, y_o) = exp[-10(x ² _o+3y ² _o)+j0.2sin(48πx ² _o)cos(32πy ² _o)], which has a 2D Gaussian function but with different profiles and chirp frequencies in each direction. Figs. 4(a), 4(b) and 4(c) and 4(d) are the original and rebuilt signals respectively. The rebuilt signal is still consistent and well agrees with the original one. The noise effect on our algorithm is also studied. Fig. 5 shows the influence of the noise to our retrieved signal when different extent of random noise is added on our signal f(x_o) = exp[-16x ² _o+j0.5sin(64πx ² _o)] which is used as an example in Fig. 3. It can be found from the Figs. 5(a) and 5(b) that the recovered signal is good for signal to noise ration (S/N) is equal to 10, but it would have large errors when the S/N is equal to 3. Thus this method has good noise resistance property.

 

Fig. 3. Numerical results of the proposed recursive algorithm. (a) The magnitude of the recovered signal (the open circles) and the magnitude of original signal f(x_o) (the solid curve); (b) The phase of the recovered signal (the open circles) and the phase of original signal f(x_o) (the solid curve).

Download Full Size | PDF

 

Fig. 4. Numerical results of the proposed recursive algorithm for a two dimensional signal f(x_o, y_o) = exp[-10(x ² _o+3y ² _o)+j0.2sin(48πx ² _o)cos(32πy ² _o)]. (a) The magnitude of the original signal (b) The phase of the original signal (c) The magnitude of the recovered signal (d) The phase of the recovered signal.

Download Full Size | PDF

 

Fig. 5. Numerical results of the proposed recursive algorithm with different extent of random noise. (a) and (b): The magnitude and phase of the original signal (the solid curve) and the recovered signal (the open circles) with S/N = 10 (c) and (d): The magnitude and phase of the original signal (the solid curve) and the recovered signal (the open circles) with S/N = 3.

Download Full Size | PDF

5. Summary

The discrete Fresnel transform (DFrT^z) and its essential properties were presented in this paper. A new concept called “complex convolution” relation was introduced. Using both the correlation and the complex convolution relations, we developed an improved recursive algorithm that can reconstruct a symmetrical signal using only one Fresnel transform intensity. The numerical results showed that this scheme is successful and both the signal’s magnitude and phase information can be rebuilt very well. The advantages of this method are that it is very efficient and only one Fresnel-intensity is required, but it can only be applied to symmetric signals or figures. Comparing it with the oversampling phasing method [10-12], the latter can be applied to any forms of images or objects and only one Fourier-transform-intensity needed which is easier to obtain the Fraunhofer image for the X-ray diffraction microscopy, although it is still an iterative method and some constraints are needed. This algorithm can be only applied to a symmetrical signal. However, any signal and its own mirror signal can be synthesized into a symmetrical signal, making this method practicable. As for the speed of our recursive algorithm, the authors have estimated carefully the amount of the arithmetic operations involved in our program. We compare it to that of another previous recursive work [9] and found that our speed should be at least twice faster. The computation time on a personal computer also showed the same trend for these two algorithms. Compared with previous works, the proposed method is advantageous in some ways; such as, no lenses are needed in the process and this method is efficient and economical to implement because only one intensity is necessary.