1. Introduction

The reconstruction of a signal from the magnitude of its Fourier transform (F.T.M.), generally referred to as the phase retrieval problem [1], arises in a variety of different contexts and applications and within such diverse fields as X-ray crystallography [2], X-ray imaging[3, 4], electron microscopy [5], astronomy, optics, and signal processing [1,2,3,6]. This problem, however, does not have a unique solution in general [1]. To find a solution for this problem, researchers have tried several methods. One approach requires additive signals. For example, Taylor proposed an off-axis holography technique by adding a point signal far off the region of the support of the desired signal [7]. Fienup considered the reconstruction of signals having latent reference points [8]. Fiddy et al. proposed a method using Eisenstein’s criteria to make a two-dimensional signal irreducible [9]. Kim and Hayes sought phase retrieval by adding an arbitrary signal to an unknown desired signal and developed several conditions and reconstruction algorithms [10]. Kim showed that if an additive point signal having large amplitude was added to the origin of a one-dimensional signal, the signal becomes a minimum-phase signal so that the signal can be uniquely defined from the magnitude of its Fourier transform. In addition, if a point signal having large enough amplitude is added to the origin of a two-dimensional signal, the added signal can be decomposed into one-dimensional minimum-phase signals that can be determined from the magnitudes of their Fourier transforms [11].

Among these, [10] explains the problem of reconstructing an unknown one- or two-dimensional signal from the magnitudes of the Fourier transform of two signals. The unknown desired signal and the other is a signal given by the addition of the unknown desired signal and a known “reference” signal. The authors presented several conditions under which the signals can be uniquely specified from the magnitudes of the two Fourier transforms and presented closed-form algorithms which enable us to reconstruct the desired signals under uniqueness conditions.

In this paper, as a continuation of the research in [10], we discuss phase retrieval when the additive signal is a point signal and the desired signal is one-dimensional and not specified uniquely from the given conditions. Here, there are two signals that satisfy the given conditions. We derive a closed-form algorithm that reconstructs two desired signals.

This paper is organized as follows. In Section 2, we discuss the theoretical background related to this work. In Section 3, we consider the derivation of the reconstruction algorithm. Finally, we present some examples that show the performance of the algorithms.

2. Theoretical background

In mathematical terms, the reconstruction of a one-dimensional signal x(n) from two sheets of FTM, i.e., the FTM of the signal x(n) and the FTM of another signal y(n) given by the addition of x(n) and a known point signal Aδ(n - n ₀), i.e.,

where A and n ₀ are the amplitude and the position of the point signal. To avoid ambiguity of translation and multiplication by a complex constant with unit magnitude, we assume that the original desired signal is a real, finite length discrete-time signal with the minimum nonzero support R(N) = [0,N-1] so that x(n) = 0 outside the support [0,N-1], x(0)≠0, and x(N -1) ≠ 0.

This problem looks similar to the off-axis holography technique and may be considered a generalization of the off-axis holography technique. The off-axis holography technique attempts to reconstruct the unknown, desired signal x(n) from the FTM of y(n) = x(n) + Aδ(n-n ₀) when the additive point ‘reference’ signal is located far off the region of support of the desired signal. More specifically, if n ₀ is greater than 2N - 1 or less than -N +1, the signal can be determined uniquely from the FTM of y(n) only (Please refer to Fig. 1). If the point signal is located at the origin, i.e., n ₀ = 0, and its amplitude is larger than the sum of the absolute values of the desired signal, the signal y(n) becomes a minimum-phase signal and can be determined uniquely from its FTM [11]. On the other hand, if the point signal is located at n ₀ = N -1 and its amplitude is larger than the sum of the absolute values of the desired signal, the signal y(n) becomes a maximum-phase signal and can be determined uniquely from its FTM, too [11]. Otherwise, the desired signal can not be determined from the FTM of y(n) only, but it can be determined uniquely with the help of the FTM of x(n) [10]. In [10], the authors showed that for a one-dimensional signal, this problem has a unique solution unless n ₀ = (N -1)/2 and developed closed-form reconstruction algorithms. It says in Theorems 1 and 2 in [10], if n ₀≠(N-1)/2 , the sequence x(n) is uniquely defined by the magnitude of the 2N-point DFT’s of x(n) and y(n).

 

Fig. 1. The original signal and an additive reference signal (Media 1).

Download Full Size | PDF

In this paper, we discuss n ₀ = (N -1)/2 . Since n ₀ should be an integer, we assume that N is an odd integer. There are only two signals that satisfy this condition. To justify this, let x̂(n) be another signal that satisfies the given conditions. From the results of Theorem 1 in [10], x̂(n) should be either x̂(n) = x(n) or x̂(n) = x(2n ₀ - n) = x(N -1- n) [10,12]. Since N is an odd integer, the support of x̂(n) = x(2n ₀ - n) = x(N -1- n) is [0, N - 1] , which is exactly the same as that of x(n) . Both x(n) and x(N -1- n) are solutions.

3. Reconstruction

In this section, we develop a closed-form reconstruction algorithm. The algorithm is basically the reconstruction of the original signal x(n) from its autocorrelation function and had been first used in [13]. The autocorrelation of a signal x(n) is defined as

and is given as

for a finite-length discrete-time signal with support [0, N -1]. Since the Fourier intensity, i.e., the square of the magnitude of the Fourier transform, ∣ X(e^jω)∣, of a signal x(n) and its autocorrelation r_x(k) = x(k)^* x(-k) form a Fourier transform pair,

With y(n) = x(n) + h(n), the autocorrelation of y(n) can be written in terms of the autocorrelation functions of x(n) and h(n) as follows;

Now we define r_yxh(k) as

Since we assume that ∣ X(e^jω) ∣ and ∣ Y(e^jω) ∣ can be measured and that h(n) is known, the values of signal r_yxh(k) are all known. Now, since we assume that the reference signal h(n) is a point signal

from (5) and (6), r_yxh(k) can be also written as

First, we begin with the reconstruction of x(n ₀), the center point. This is easily determined in (8) for k = 0

Second, we consider the reconstruction of the two boundary values x(0) and x(N - 1). We can not employ a linear equation to determine this. Instead, we derive a quadratic equation with the summation and the product of the two values. The summation of these two values can be derived from (8) by inserting k = n ₀ = (N - 1)/2, i.e.,

The product of the two values can be derived from Eq. (3), the definition of the autocorrelation of x(n), i.e.,

Using (10) and (11), we can form a quadratic equation

Suppose that the two roots of (12) are t ₁ and t ₂. Then, the boundary values of the two solution signals are either x(0) = t ₁ and x(N - 1) = t ₂ or x(0) = t ₂ and x(N - 1) = t ₁.

For each of these, we can determine the remaining values. From (3) and (8), we can formulate the following linear matrix equation

for n = 1,2, …, n ₀ -1. Now we assume that

From Eq. (13), the values of x(n) and x(N - 1 - n) can be determined simply if the matrix

is invertible, i.e., x(0) ≠ x(N -1). Even if this matrix is not invertible, it is possible for the remaining values to be determined.

3.1 P that is Invertible

For P to be invertible, the determinant of the matrix should be nonzero or, equivalently, x(0) ≠ x(N - 1). In this case, the values of x(n) for n = 1, 2, …, n ₀ - 1 and N - 2, N - 3, …, n ₀ +1 can be determined by solving Eq. (13), i.e.,

This completes the reconstruction of a solution signal. To reconstruct the other solution signal, follow the same procedure described above after swapping the values of x(0) and x(N - 1) or simply time-reverse the signal, or x(N -1 - n).

Now we show an example that shows how the algorithm works.

Example 2

As in Example 1, we assume that x(n) =[1 2 3 4 5] and the point reference signal h(n) =[0 0 1 0 0]. The magnitudes of the Fourier transform for x(n) and y(n) are given as ∣ X(e^jω) ∣ =[15.0 10.9 4.3 3.7 2.6 3.0 2.6 3.7 4.3 10.9] and ∣Y(e^jω) ∣ =[16.0 11.8 4.4 2.7 2.8 4.0 2.8 2.7 4.4 11.8]. From these, we get the autocorrelations of x(n) and y(n) as r_x(k)=[5 14 26 40 55 40 26 14 5] and r_y(k)= [5 14 32 46 62 46 32 14 5], for k = -4, -3 , …, 3, 4. Since the autocorrelation r_h(k) is given as r_h(k) = δ(k), from the information given, r_yxh(k) is given as r_yxh(k)= [ 0 0 6 6 6 6 6 0 0].

Now we can begin the reconstruction. Since N = 5 and n ₀ = 2, the center point x(2) can be obtained from (10)

The summation and the product of the two boundary values are r_yxh (2) = x(0) + x(4) = 6 and r_x(4) = x(0)x(4) = 5. Using these values and Eq. (12), the quadratic equation with roots of x(0) and x(4) is given by

From this, we get either x(0) = 1 and x(4) = 5 or x(0) = 5 and x(4) = 1. We discuss the first case where x(0) = 1 and x(4) = 5 . Since x(0) ≠ x(4), the remaining values of x(1) and x(4) can be determined from (15).

Thus, we get one solution x(n) =[1 2 3 4 5].

If we follow the same procedure with x(0) = 5 and x(4) = 1 or simply reverse the solution, we can get the other solution x(n) = [5 4 3 2 1].

3.2 P that is not invertible

If P is not invertible, the determinant of the matrix is zero and x(0) = x(N - 1) . In this case, the algorithm is similar to that of the case where P is invertible. However, this is somewhat complicated, so we present some key equations.

Since we have determined x(0) and x(N - 1), the next values to be determined are x(1) and x(N - 2) . From Eq (8), we can get the sum of x(1) and x(N - 2), i.e.,

From Eq. (3) for k = N -3 , we can get the product of x(1) and x(N - 2) since

and

As we have done in (12), we can construct a quadratic equation for solving x(1) and x(N - 2) as

and if we solve this equation, we get x(1) and x(N - 2).

Once the values of x(1) and x(N - 2) are determined, we can specify the values of x(n) for n = 2,3,…,(N -1)/2-1 and n = (N-1)/2 + 1,…,N - 3 similarly as follows

As we have done before, if x(1) ≠ x(N - 2), Eq. (18) has a unique solution and all the values of x(n) for n = 2,3,…,(A'-1)/2-1 and n = (N - 1)/2 + 1,(N - 1)/2 + 2,…,N - 3 can be determined. If x(1) = x(N - 2), we can do similar work to what is presented here.

If p_n is not invertible for n =0, 1, …, N/2-1, then this means that the signal is symmetric. Now we consider an example to see the performance of the algorithm above.

Example 3

As an example, we assume the unknown original signal is given as x(n) = [1, 2, 3, 4, 5, 6, 7, 8, 1] and the known additive point reference signal is h(n) = [0, 0, 0, 0, 3, 0, 0, 0, 0]. The signal y(n) is given as y(n) = x(n) + h(n) =[1, 2, 3, 4, 8, 6, 7, 8, 1]. In this example, N and n ₀ are N = 9 and n ₀ = 4, respectively.

The values of the FTM of x(n) and y(n), which are measurable quantities, are given as ∣ X(e^jω) ∣ =[3.59, 3.36, 4.36, 4.61, 6.40, 7.94, 12.85, 28.68, 37.00, 28.68, 12.85, 7.94, 6.40, 4.61, 4.36, 3.36, 3.59, 3.00] and ∣Y(e^jω)∣= [3.68, 5.79, 7.21, 5.72, 3.64, 6.00, 14.81, 31.56, 40.00, 31.56, 14.81, 6.00, 3.64, 5.72, 7.21, 5.79, 3.68, 0]. From these, we can calculate the autocorrelation functions of x(n) and y(n), i.e., r_x(n) = [1, 10, 26, 48, 75, 106, 140, 176, 205, 176, 140, 106, 75, 48, 26, 10, 1] and r_y(n) =[1, 10, 26, 48, 81, 136, 170, 206, 244, 206, 170, 136, 81, 48, 26, 10, 1] for n = -8, -7, ⋯, 7, 8. Since the autocorrelation r_h(n) is given as r_h(n) = 9δ(n) from the given information, the function r_yxh(n) is given as r_yxh(n) =[0, 0, 0, 0, 6, 30, 30, 30, 30, 30, 30, 30, 6, 0, 0, 0, 0] for k = -8, - 7, ⋯, 7, 8.

From Eq. (10), x(n ₀) is calculated as follows:

To get the two boundary values, the quadratic equation of Eq. (12) is derived as

of which two roots are x(0) = x(8) =1. Since in this case, the matrix P in Eq. (13) is not invertible, we are not able to determine the remaining values using Eq. (15). Instead, we determine the next outer most values x(1) and x(7).

From Eq. (18) and (19), the sum and the product of x(1) and x(7) are given as x(1) + x(7) = r_yxh(3)/3 = 10 and x(1)x(7) = r_x (6) - x(0)r_yxh (2)/3 = 16, the quadratic equation (19) becomes

If we solve this equation, we get the set of solutions as either x(1) =1 and x(7) = 8 or x(1) = 8 and x(7) =1. For the first case, the matrix equation of Eq. (20) is given as

By solving this equation for n = 3, 4, we get x(2), x(3), x(6) , and x(5), which makes the entire solution

For the second case, we get the other solution signal either by deriving and solving the equations (19)–(20) or by simply reversing the solution above, i.e.,

In Fig. 2, we show the block diagram of the reconstruction algorithm.

 

Fig. 2. Block Diagram of the Reconstruction Algorithm (Media 2).

Download Full Size | PDF

4. Noise analysis

Since the reconstruction equations (15) and (20) are composed of recursive algorithms, this algorithm is expected to be sensitive to various noises including additive Gaussian noise and round-off errors occurring during the calculation. In this section we consider the effect of various noises and parameters. All the calculations are done in Matlab and in double precision. The mean square error is defined as

4.1 Effect of the length of the original signal

In this experiment, we changed the length of the original signal from 7 to 127. The MSE is defined as Eq. 21 and shown in Fig. 3 (a). As the length of the original signal increases, the number of computation also grows and computational noises due to truncation error also grow. As is expected the errors are growing as the size of the support increases. However, the maximum error is less than 1 and thus the size of the support does not affect much.

4.2 Effect of the amplitude of the point function.

The amplitude of the point function also can affect the error. In this experiment, we changed the value of amplitude, i.e., A in Eq. (1), from 1 to 10,000 and assumed that the support of the original signal is given as R[0,127]. The results are shown in Fig 3(b). As we can see in this graph, the error decreases as the amplitude increases but does not decrease further once the error reaches a certain level.

4.3 Effect of the additive Gaussian noises

To see the effect of additive Gaussian noises, we added Gaussian noises to the measured Fourier transform magnitudes. We assumed that the length of the original signal is 127 and the standard deviation varies 0 to .1 by the incremental step of 0.005. As we can see in the Fig. 3 (c), the mean square error grows fast once the noises are added and however the error does not grow much as the standard deviations of the noises increase. The left most is the case when there is no noise. To see the characteristics of the errors, in Fig 3(d), we showed the original signal (blue solid line), reconstructed signal (red solid line), and the difference signal between the two signals (blue solid line with points at the bottom) when the standard deviation of the additive Gaussian noise is 0.1. The difference values, i.e., the error values, become bigger near the center of the support of the signal. This is somewhat expected results since the algorithm is composed of recursive equations and the reconstruction is done from the end to the center of the support. We can see in this graph that although the overall MSE is big, the error is concentrated at near the center of the signal and there is not much error near the end points.

 

Fig. 3. Effect of noises and other parameters on the mean squared error: (a) (Media 3) MSE versus the length of the original signal, (b) (Media 4) MSE versus the amplitude of the additive point signal, (c) (Media 5) MSE error versus standard deviation of the additive Gaussian noises, and (d) (Media 6) An example of the original (blue solid line), reconstructed (red solid line), and the error signal (blue solid line with points at the bottom).

Download Full Size | PDF

5. Summary

In this paper, we considered the problem of reconstructing an unknown one-dimensional signal from its FTM and the FTM of another signal that is related to the unknown signal through the addition of a known point ‘reference’ signal at the center of the non-zero support of the desired signal. We show that there are two signals that can be specified from the two FTMs. We presented closed-form recursive algorithms that may be used to reconstruct these two solution signals. We hope that the results of this work may be helpful in solving real-world problems related to optics, astronomy, optical communications, and X-ray crystallography.