Abstract
Inspired by the implementation of the fractional Fourier transform (FRFT) and its applications in optics, we address the problem of reconstructing a signal from its several FRFT magnitudes (or intensities). The matrix completion method is adopted here. Through numerical tests, the matrix completion method is proven effective in both noisy and noise-free situations. We also compare our method with the Gerchberg-Saxton (GS) algorithm based on FRFT. Numerical tests show that the matrix completion method gains a certain advantage in recovering uniqueness and convergence over the GS algorithm in the noise-free case. Furthermore, in terms of noisy signals, the matrix completion method performs robustly and adding more measurements can generally increase accuracy of recovered signals.
1. Introduction
Phase retrieval has a long history since the first application to the optics and X-ray crystallography [1]. And in the last few decades phase retrieval has arisen in various fields including diffraction imaging [2], image encryption [3], and quantum mechanics [4]. In the common case, it is not enough to recover the signal only with the magnitudes of the original signal and its Fourier transform. To overcome nonuniqueness or improve the probability of successful recovery, two approaches are usually adopted. The first is to exploit extra prior information, like the support of the signal. The second approach is to increase the number of measurements, for example, oversampling and using the short-time Fourier transform and multiple linear canonical transforms. For more detail, we refer readers directly to [5–8]. In this paper, we address the problem of recovering a signal from its multiple FRFTs in different orders, which can be categorized into the second approach. Multiple FRFT magnitudes can provide more information of the signal, which may lead to the correct recovery without extra constraints of the original signal. This is a notable advantage because one of most stringent limitations for most algorithms is the need for isolated objects (the support constraint) [9].
Furthermore, phase retrieval from magnitudes of FRFT is a problem of practical interest and has been studied extensively in many areas, such as image encryption and X-ray tomography [10–12]. A direct application is in ABCD optical systems, which consists of thin lenses and sections of free space. In such systems, the transverse amplitude of light of different planes can be related through FRFT. In other words, the method proposed here can be used to recover complex light field from transverse amplitude of light at several planes. Phase retrieval also has a profound application in image encryption [10]. One basic step in image encryption based on FRFT is to recover the original signal accurately from FRFT intensities of encrypted images. And, in this process, using a stable and accurate recovering methodology is of crucial importance.
On the other hand, phase retrieval based on FRFT is a direct generalization of the common phase retrieval based on Fourier transform. In quantum physics, Pauli [13] first put forward the question whether it is sufficient to uniquely determine the wave function with the knowledge of position distribution and momentum distribution . can be understood as the Fourier transform of . For a set of magnitude measurements, phase retrieval is unique if any two solutions and can be related by , where the norm of the scalar is 1. Unfortunately, without additional information or constraints, and do not uniquely determine up to a global phase. Then it is natural to consider recovering a signal from multiple FRFT magnitudes.
As for how many FRFTs for phase retrieval suffice for unique recovery, there are some valuable results. In the continuous case, Millane explains shortly in [14] why conventional phase retrieval can cause ambiguity. Carmeli et al. elegantly proves that three FRFT magnitudes, regardless of angles chosen, are not enough to guarantee the uniqueness of phase retrieval of arbitrary signals in continuous case. Whether four FRFT magnitudes suffice or not is still an open question [15]. As for the discrete case, it is proven that observables of a signal of length can ensure recovering uniqueness with a great probability [16]. Therefore, we consider using no less than four discrete FRFTs to ensure the uniqueness of recovery. It should be noted that how to choose FRFTs to make recovery unique has not been solved yet. However, phase retrieval from more than two FRFTs has not been paid much attention to, except for a few of articles [17, 18]. Currently, most studies deal with magnitudes of two FRFTs, which need some additional assumptions about the signal, like support of the signal and nonnegativity. Without these assumptions algorithms dealing with two FRFTs cannot ensure recovering correctness in many cases [16, 19]. In this paper, at least 4 FRFTs are chosen to reconstruct the original signal to ensure probability of successful recovery without further assumption upon signals.
Then it comes to the recovering algorithm. Many algorithms have been proposed to solve the problem of phase retrieval like the iterative algorithm [20, 21] and the Gerchberg-Saxton (GS) algorithm [22]. The GS algorithm possesses a dominant role and a large number of algorithms are developed on the basis of the GS algorithm. Phase retrieval from magnitudes of multiple FRFTs using the GS algorithm has been proposed and comprehensively studied by [17]. In this paper, the matrix completion method [19] is adopted to recover the signal from magnitudes of FRFTs, which has not been researched before. We transform the problem into a matrix completion problem and then solve it by the convex optimization algorithm.
The matrix completion method is tested both in the noise-free case and in the noisy case. It should be noted that no extra prior assumption like support constraint is needed for the signal. And numerical results verify the feasibility of this method. We also compare the recovering convergency of the matrix completion method and the GS method. Numerical experiments show that the matrix completion method gains certain advantages in stability and convergency.
2. Matrix Completion Method for FRFT Phase Retrieval
The FRFT is usually defined as [23] is a kernel function defined aswhere is some integer, , and is the Dirac function. FRFT has the following important properties. (i) When , FRFT turns into the ordinary Fourier transform. (ii) FRFT is additive in index; namely, . (iii) FRFT is a unitary operator.
Since we always deal with the digital data in practice, it is necessary to introduce the corresponding discrete fractional Fourier transform (DFRFT). There are many types of DFRFT [24–26]. The GS method in [17] is based on eigendecomposition-type DFRFT. Therefore, we also use this type of DFRFT proposed in [26] for comparability in numerical experiments. The DFRFT matrix is generated through calculating the power of the discrete Fourier matrix, which involves the discrete Hermite function. However, the discrete Hermite function cannot be expressed in analytical form. Therefore, we calculate the DFRFT matrix numerically according to [26]. Besides, the eigendecomposition-type DFRFT can preserve all important properties of continues FRFT and approximates FRFT well, despite the lack of analytical expression.
The matrix completion method is stated as follows. Suppose we do times FRFT separately upon the signal , and the square of magnitudes of transform is accordingly . Then the knowledge of square FRFT magnitudes can be written in the following forms of quadratic constraints:where , are DFRFT matrix of different orders; the th row is denoted as . is the th element of .
Then we transform this problem into a semidefinite program by lifting up the quadratic constraints to measurements about rank-one matrix: . Specifically,where . Here we use operator to define the map from positive semidefinite matrix into . For simplicity, we defineThen the phase retrieval problem can be rephrased asHowever, this problem is proven to be NP hard. In fact, (6) just rephrases the original problem and it does not make the original problem easier to solve. To make it trackable, a relaxation of (6) is adopted here. The relaxed form is a convex optimization problem as follows:Equation (7) can be solved by various numerical methods such as Nesterov’s method [27]. Currently, there is still no proof of the equivalence of (7) and (6). However, empirical numerical results indicate the feasibility of (7). Assuming is the solution of (7), we can then take as estimation of , where is the best rank-one approximation of . Practically, (7) is transformed into the following convex optimization problem:where is a real positive number. The corresponding matrix completion algorithm is presented in Algorithm 1.
Algorithm 1. The matrix completion algorithm for recovering the signal from multiple FRFT magnitudes is as follows:
Input. FRFT square magnitudes .
Output. Estimate of original signal (i)Obtain by solving (ii)Return , where is the best rank-one approximation of
3. Numerical Experiments
In this section, we illustrate the effectiveness of the matrix completion. We test our matrix completion method in both the noise-free case and noisy case, and a comparison between the matrix completion method and the GS method is also made.
3.1. Numerical Solver and Setup
The core part of Algorithm 1 is solving a convex optimization problem (9). There are many numerical solvers for solving this problem and we choose TEFOCS [28] in this paper. TEFOCS adopts a projected gradient method which starts from an initial guess . We initialize by setting as a zero matrix. The TEFOCS terminates if two successive iterate results, and , satisfy where denotes the Frobenius-norm of matrix. And the threshold value is taken as in numerical experiments. In terms of the positive scalar , one can refine the algorithm’s performance through selecting the value of to make and have close orders of magnitude [27]. Therefore, in noise-free case, is simply taken as a small value close to and, in noisy case, is larger.
To judge the error between the recovered signal and the original signal , we use the relative mean-square error (RMSE) which is defined as In practice, we could calculate RMSE by . In addition, all the numerical experiments in this section are based on double-precision floating-point numbers. However, practical platforms such as FPGA and VLSI do not support floating-point numbers [29–31]. Therefore, improving accuracy in engineering applications might require more investigations in future.
3.2. Noise-Free Situation
For noise-free signals, we test our method on two different signals with the same length and take . The first signal is a smooth complex signal which consists of several sinusoids. The second signal is a complex signal with each element independently obtained from uniform distribution (for each entry , and obey independent uniform distribution over ). The orders of discrete FRFTs are . The recovering results are presented in Figure 1. We can see that the original signal and the recovered signal are indistinguishable to eyes.
(a) Smooth signal and its recovery (real part)
(b) Random signal and its recovery (real part)
3.3. Noisy Situation
We also evaluate the matrix completion method in terms of noisy signals. The parameter is simply taken as 2 empirically. However, it is important to emphasize that in practice one would better find an ideal value via a strategy like generalized cross validation (GCV) [19]. We add several Gaussian noises under different intensities to the signal and reconstruct it from different numbers of magnitudes of FRFT. The test signals are random signals obeying uniform distribution over . In each SNR level we repeat the experiment 15 times. It can be seen in Figure 2 that RMSE decreases almost linearly as SNR increases, which shows that the matrix completion method is robust with noise. And, as for the recovering accuracy, Figure 2 shows that RMSE of the method using less FRFTs is obviously larger than that of method using more FRFTs. Therefore, one can improve accuracy of recovering results by adding measurements. Nevertheless, we can see that the RMSE lines of 6-FRFT and 7-FRFT method are very close, which indicates that adding measurements may not improve the accuracy markedly when the number of measurements is large enough. However, it is worth noting that the accuracy can be improved if we utilize the noise’s property.
3.4. Comparison with the GS Algorithm
Lastly, we compare our method with the GS method [17]. Both these two methods do not need prior assumption of signals.
Algorithm 2. The GS algorithm for recovering the signal from multiple FRFT magnitudes is as follows: Input: , max iter, FRFT magnitudes and initial guess . Output: Estimate of original signal for to max iter do for do , for end for if then Break end if end for return
The GS algorithm to recover from multiple FRFT magnitudes is illustrated in Algorithm 2. Suppose we have the multiple FRFT magnitudes of : . We initialize . The stopping criterion is the fact that is less than a given threshold value, , which is taken as in this paper, or iteration times are more than max iter, which is given as . The latter condition is necessary in case the GS algorithm does not converge.
Although the iterative procedure of the GS algorithm is simple and efficient, the convergence of it cannot be guaranteed more probably. We test the GS algorithm upon the signal , , and the sampling rate is 50. The chosen orders of DFRFTs are, respectively, . Figure 3(a) demonstrates how RMSE of the GS algorithm varies at each iteration step. When , the signal recovered by the GS algorithm does not converge to the original signal and the RMSE stays at a certain level after some iteration steps. Then we add measurements and the RMSE can diminish when . However, convergence does not remain when is 6 or 7; the RMSE is still quite considerable. By contrast, the matrix completion method is also tested and the RMSE curves can all steadily vanish to zero as shown in Figure 3(b). In terms of the computational complexity, the matrix completion method is considerably less efficient than the GS method, since it needs to deal with a relatively large-scale convex optimization problem. However, the matrix completion method has an important strength in convergence to the original signal.
(a) RMSE after each iteration (GS)
(b) RMSE after each iteration (matrix completion)
To judge the convergence more carefully, we implement these two algorithms 200 times repeatedly upon random signals in length of 32: , where and are independent uniform distribution variables over . If RMSE is larger than 10%, we declare that the reconstruction fails. The failure rates of these two algorithms are illustrated in Table 1. It is interesting to see that the GS failure rate decreases as increases. However, the matrix completion method performs better with a lower error frequency.
4. Conclusions
In this paper we consider solving the FRFT phase retrieval problem with the matrix completion method which utilizes magnitudes of multiple FRFTs. In numerical tests we can see that this method works well for both noisy signals and noise-free signals. Compared with the GS method, the matrix completion method has a larger probability in recovering the original signal, especially when there are less FRFT magnitudes. This method also performs robustly in the presence of noise. Besides, one can improve recovering accuracy as well as probability through adding more FRFT magnitudes. However, this conclusion still needs to be confirmed by further rigid theoretical investigation. In future we will find a practical strategy to improve the efficiency of the matrix method.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This project is supported by the National Natural Science Foundation of China (Grant no. 61571008).