Signal reconstruction algorithm based on a single intensity in the Fresnel domain

A novel algorithm that can reconstruct a symmetrical signal (both the amplitude and the phase information) with only a single Fresnel transform intensity is proposed. A new complex-convolution method is introduced, which is needed in the algorithm. The essential properties of the discrete Fresnel transform are presented as well. Numerical results show that this method can successfully rebuild the signal from one signal intensity, which is more advantageous in speed and efficiency than the conventional method that requires two intensities to accomplish this task. ©2007 Optical Society of America OCIS codes: (070.4560) Optical data processing, (100.2000) Digital image processing, (200.4740) Optical processing. References and links 1. H. A. Ferwerda, “The phase reconstruction problem for wave amplitudes and coherence functions,” in Inverse Source Problems in Optics, H. P. Baltes, ed. (Springer-Verlag, Berlin, 1978) 2. 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) 3. R. W. Gerchberg and W. O. Saxton, “Phase determination for image and diffraction plane pictures in the electron microscope,” Optik 34, 275-284 (1971) 4. P. Van Toon and H. A. Ferwerda, “On the problem of phase retrieval in electron microscopy from image and diffraction pattern. IV. Checking of algorithm by means of simulated objects ,” Optik 47, 123-134 (1977) 5. Z. Zalevsky and R. G. Dorsch, “Gerchberg–Saxton algorithm applied in the fractional Fourier or the Fresnel domain ,” Opt. Lett. 21, 842-844 (1996) 6. W. J. Dallas, “Digital computation of image complex amplitude from image and diffraction intensity: an alternative to holography,” Optik 44, 45-59 (1975) 7. W. Kim and M. H. Hayes, “Phase retrieval using two Fourier-transform intensities,” J. Opt. Soc. Am. A 7, 441-449 (1990) 8. N. Nakajima, “Phase retrieval from two intensity measurements using the Fourier series expansion,” J. Opt. Soc. Am. A 4, 154-158 (1987) 9. W. X. Cong, N. X. Chen, and B. Y. Gu, “Phase retrieval in the Fresnel transform system: a recursive algorithm,” J. Opt. Soc. Am. A 16, 1827-1830 (1999) 10. C. Song, R. J. Damien, Y. Nishino, Y. Kohmura, T. Ishikawa, C. C. Chen, T. K. Lee, and J. Miao, “Phase retrieval from exactly oversampled diffraction intensity through deconvolution,” Phys. Rev. B. 75, 012102 (2007) 11. M. A. Pfeifer, G. J. Williams, I. A. Vartanyants, R. Harder, and I. K. Robinson, “Three-dimensional mapping of a deformation field inside a nanocrystal,” Nature 442, 63-66 (2006) 12. J. Miao, P. Charalambous, J. Kirz, and D. Sayre, “Extending the methodology of X-ray crystallography to allow imaging of micrometer-sized non-crystalline specimens.” Nature, 400, 342-345 (1999) 13. J. Miao, D. Sayre, and N. H. Chapman, “Phase retrieval from the magnitude of the Fourier transform of nonperiodic objects,” J. Opt. Soc. Am. A. A15 1662-1669 (1998) 14. Y. M. Bruck and L. G. Sodin, ‘‘On the ambiguity of the image reconstruction problem,’’ Opt. Commun. 30, 304–308 (1979) 15. M. H. Hayes, ‘‘The reconstruction of a multidimensional sequence from the phase or magnitude of its Fourier transform,’’ IEEE Trans. Acoust., Speech, Signal Process. ASSP-30, 140–154 (1982) (C) 2007 OSA 2 April 2007 / Vol. 15, No. 7 / OPTICS EXPRESS 3766


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][2][3][4][5][6][7][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][11][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.

Theory
The Fresnel transform ( z FrT ) of a function ( ) where z is the propagation distance, λ is the wavelength, o x is the variable of the coordinate of the input plane, and p x is the variable of the coordinate of the output plane.As illustrated in Fig. 1, the optical field in the output plane p p x y is ( )

∑
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 .This must satisfy the Nyquist sample theory that o x δ 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][14][15]. The where ( ) 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  5) and ( 6), the correlation property can be rewritten as where ( 2) , and let ( ) R k be the entire right-hand side of Eq. ( 7) as If the autocorrelation lag is set as 1 k N = − , from Eqs. ( 7) and ( 8), we have Next, for the lag .
Similarly, we derive the relation, for the lags 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 = − , 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 − ) 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 ( ( ) 12), and let ) (k R′ be the entire right-hand side of Eq. ( 12) as We have the following equation as Next, for .
For the lags k N m = − + with 2 m ≥ , 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 1 k N = − + , and so on; and the symbol B(i,j) is used to identify the indicated term.Again, only B(1,1) ) in the first row can be obtained directly by Eq.( 14), but these product factors are very important for our aim as explained later.

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) ) are calculated from Eqs. ( 9) and ( 14).It is found that B(2,1) 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 * ( 2 1) ( 2) With the help of the symmetrical condition ( 2 ) can be obtained when A(2,2) ) is substituted in Eq. (10).Once this is found, the factors B(3,1) ) 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) ) and B(4,3) Therefore (0) f 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 (0) f is known and refer to the first term and last term from bottom to head in part A. We can divide * ( 2 ).And we can also divide * (0) ( 2 1) ) f N − .For the same reason, the value of (1)  f can be found when * ( 2 ).This process is repeated and all the sequence [ ] can be calculated.The signal is thus reconstructed.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 ( ) 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 , 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 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.S e t in i t i a l I 0

Summary
The discrete Fresnel transform ( z DFrT ) 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][11][12], the latter can be applied to any forms of images or objects and only one Fourier-transformintensity 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.

2 .
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 3 k N = − to solve the factor A(3the 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.This recursive process continues until the value ( shown in part A triangle) and every product term in partA is found.Finally we can take the product of the two known terms * part A triangle respectively) and use Eq.(9) to get.

Fig. 3 .
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 ( ) o f x (the solid curve); (b) The phase of the recovered signal (the open circles) and the phase of original signal ( ) o f x (the solid curve).

Fig. 5 .
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. inverse N − .With the help of Eqs. ( is on the second row and the first term.