CC BY
37Electronic Journal «Technical acoustics» http://webcenter.ru/~eeaa/ejta/
2 (2002) 15.1-15.6 Z. Chama1, M.E. Louhibi1, M.F. Belbachir2, A. Djebbari1
1 Telecommunication and Digital Signal Processing Laboratory Electronic Department, University of Djillali Liabes in Sidi Bel Abbes BP 89 22000, Algeria
2 Signal and Systems Laboratory University USTO, Oran, Algeria
Two dimensional signal reconstruction from its spectral phase
Received 14.08.2002, published 19.12.2002
In this paper, a new algorithm of image reconstruction from the phase of its Fourier transform is presented. It uses the eigenfilter method to approximate a given phase response in the least square sense. This method has proved its efficiency as the iterative methods. An example will be given to show the importance of this algorithm and the high image reconstruction quality.
1. INTRODUCTION
The Fourier transform phase or magnitude information alone is not, in general, to uniquely specify a signal, the ability to reconstruct a signal from only phase or magnitude information would be useful in number of important applications. For example, in many problems which arise in X-ray crystallography, electron microscopy, coherence theory, and optics. In many applications, either the spectral magnitude or phase of a signal may be severely distorted component. For example, in the class of problems referred to as blind deconvolution, a desired signal is to be recovered from an observation, which is the convolution of the desired signal with some unknown signal [1]. Since little is usually known about either the desired signal or the distorting signal, deconvolution of the two signals is usually a very difficult problem. However, in the special case in which the distorting signal is known to have a phase, which is identically zero, the phase of the signal is undistorted. Such a situation occurs, at least approximately, in long-term exposure to atmospheric turbulence or when images are blurred by severely defocused lenses with circular aperture stops [2]. In this case, except for phase reversals, the phase of the desired signal and, therefore, it is of interest to consider signal reconstruction from phase information alone. Several iterative or recursive methods have been proposed for this purpose. Based on the Fourier transform [3, 4, 5, 6-9], the short time Fourier transform [10] or the bispectrum [11, 12].
An algorithm under which an image is fully specified within a scale factor from the phase of its Fourier transform is developed in this paper.
It has been effectively shown that a discrete time signal may be fully specified by the phase of its Fourier transform under some specific conditions (not simply minimum or maximum phase): the zeros of its z-transform must not be located on the unit circle nor occur in complex reciprocal pairs [4]. We present the proposed reconstruction algorithm in section 2, and give numerical examples in section 3, and section 4 concludes with a summary.
2. DESCRIPTION OF THE ALGORITHM
In this section, the problem of reconstructing a multidimensional sequence from the phase of its Fourier transform is addressed. In order to ensure that there is a unique solution to the reconstruction problem, it will be assumed that x(n) is a m-dimensional sequence with support R(N) which has a z-transform with no symmetric factors. One algorithm for reconstructing x(n) from the phase of its Fourier transform involves finding the solution to a set of linear equation. Specifically, from the definition of )x (a), it may be shown, as in [4], that x(n) satisfies the equation
^ x(n) sin [)x (a) + noi] = -a0 sin )x (a),
n^Q
(1)
n K
provided that )x (a) ^± —. For the case in which )x (a) = ± —, then x(n) satisfies the
equation
^ x(n) = cos (na) = -aQ.
n^Q
(2)
Substituting the values of the phase of the M0 -point discrete Fourier transform into (1) leads to a set of M0 linear equations in N0 unknowns where M0 = M1 xM2 x....xMm and N0 = N1 xN2 x....xNm. Due to the symmetry of the phase, half of these equations are redundant and may be eliminated. Arranging the elements of x(n) into a vector of length N0 , x(n)={x(1X x(2)........., (N0)} } is a (N0 xM0) matrix, and b is a vector of length M0.
' sin[) («)+«] sin[)2 («)+«] • sin[)MQ (a)+a] " ' sin)(«) ^
A= sin[) (a)+ 2a] sin[)2 (a)+ 2a] • sin[)MQ (a)+ 2a] , b= sin)2 (a)
sin[) (a)+ NQa] V sin[)2 (a)+ NQa] • • sin[)MQ (a)+ Noa] Q y sin)M o(a) V Q /
The linear equations in (1), when augmented with the equation x(0)=a0, may be written
as
Ax = -aQ b.
(3)
The solution of (3) corresponds to the reconstructed sequence Vx, which has a M0 -point DFT, which phase tangent equals that x(n). Therefore, the matrix S = AT A, is non-singular and the desired sequence may be reconstructed from
Vx = -a„ ((A)ta)"' (A)-1 b. (4)
Although this algorithm has been used successfully in reconstructing 2-D sequences of moderate size (N0 < 256), its application in practice is limited by the computational difficulties encountered in solving linear equations when the number of unknowns becomes
large. Therefore, it’s of interest to consider alternative solutions to the reconstruction problem.
The problem may be viewed as the design of all-zero filter using a least square (LS) criterion, by the eigenfilter method, given a prescribed phase response [13]. This technique was used to reconstruct a signal from its Fourier transform phase [14]. A review of the algorithm is presented in one-dimensional form.
The z-transform of a discrete time sequence x(n) , n=0, 1... N-1, is
X(z) = £x(k)z.
(5)
k=0
Let z=e]<°, by introducing the column vectors x, C and S, whose elements are xk(a) = x(k), Ck(a) = cos(ka), Sk(a) = sin(ka), where k = 1, 2,...,N.
The Fourier phase of this sequence can be expressed as in Eq. (6).
0x (a) = arg( (e](0^)= - arctan
IX —I
^ x(k) sin(ka)
k=0
IX —i
^ x (k) cos(ka)
k=0
=- arctan
xTST
xTCT
(6)
where the superscript T means transposition.
Consider a prescribed phase dpre (a). The least square method consists to minimize the following LS error function (7).
ji
Els =jw (a)| A0(a)|2 da,
(7)
where AO (a) =0pre (a) -0X (a) is the phase error shaping function, and W(a) is the
weighting function.
Using Eq. (6), the error function (7) may be rewritten as in Eq. (8).
E
LS
arctan
( xTSP (a) A xTCP (a)
/-j
da
(8)
where SP (a) and CP (a) are the column vectors with elements
SPK (a) = sin(Opre (a) + ka), Cpk (a) = cos (Opre (a) - ka), where k = 1, 2, ..., N.
Since the error function must be kept as small as possible, a straightforward approximation is to replace the arctangent function by its argument (x ^ 0, arctan (x) ~x), as in Eq. (8). The approximate LS error is then given by Eq. (9).
E r xTSp(a)SpT(a)x da.
LS 0 x CP (a)Cp (a)x
(9)
Despite the simplification, a tricky non-linear optimisation remains to be performed. Let us, however, consider the ideal case of perfect reconstruction, the integrand in Eq. (9)
2
nullifies at all frequencies a , considering only the optimum of numerator, which is also nullifies. A new simpler quadratic error function may thus be defined as in Eq. (10).
Eqls = x Px, (10)
where P is the real-valued, symmetric and positive definite matrix of order N by N defined by Eq. (11):
ji
P = J SP (rn)STP (rn)drn. (11)
0
Now, the problem may be expressed as follows. Find the vector x such that the quadratic form Eqls is minimised, subject to the constraint xTx = 1. The solution vector xx is clearly
the eigenvector corresponding to the minimum eiguenvalue of matrix P, according to Rayleigh’s principle [15]. The validity of using expression in Eq. (10) in place of Eq. (9) must be tested in practical cases.
3. NUMERICAL EXAMPLES
The phase spectrum may be obtained via the discrete Fourier transform (DFT) algorithm. In this case, the length M0 must be greater or equal twice the signal length N0, since we
keep only the phases for the positive frequencies (a e [0, n]), i.e. half the number of the DFT points.
The eigenfilter algorithm described by Eqs. (10) and (11) is then applied to the phase responses (considered as the described )pre). For each value of M0, the reconstructed
sequence xx is determined. The value of the factor /5 is obtained using one of the following alternative definitions:
N0 /N0
5 = X x(k)/ X xs (k) (mean criterion),
k=1 / k=1
N0 /N„
15 = X |x(k)|/ X |xs (k)| (L1 norm criterion), (12)
/5 = x(k)/xs (k) for some index value k (direct comparison criterion).
In the case of 2-D signal (m=2, example image) theses criterions are not sufficient to achieve a perfect reconstruction, therefore, it is preferably to use the correlation criterion to determine the quality of reconstruction by considering the range [0,1], (1: corresponds to a perfect reconstruction, 0: bad reconstruction)
An example in Figures 1 and 2 is presented. In Fig. 1(a), a binary image of 20x32 pixels is shown, which is arranged as 1-D vector of length 640 elements. In our algorithm it is supposed that the available information is the phase, then the image will be reconstructed from its phase using a 1280-point DFT. A small size image has been chosen to avoid an expensive computational complexity. The image in Fig. 1(b) shows the phase-only representation of this image, obtained by setting the DFT magnitude equal to a constant. Using the eigenfilter method the reconstructed image is mentioned in Fig. 2 (image has been scaled for display).
(a)
(b)
Fig. 1. Original image and its phase-only synthesis (a) original image, (b) phase-only image formed by combining the phase of the Fourier transform of image (a) with a constant magnitude
E denotes the total error defined by Eq. (13) between the estimate and original sequence xe (n) and x(n) respectively, the image in this example is considered as a one-dimensional
vector by transforming the matrix image in one-dimensional vector. E is found to be equal to 10-20, and the 2-D correlation coefficient, which is equal to 1, that it is significant to a perfect reconstruction and a high reconstructed image quality.
iv —1
E = X[(n) - Xe (n)] •
(13)
Fig. 2. Phase-only image reconstruction. Image reconstructed from phase-only image (a) original image, (b) reconstructed image
4. CONCLUSION
A new algorithm for image reconstruction from the phase of its Fourier transform has been elaborated, which does not rely on a minimum or maximum phase property. The idea is to design an all-pass filter with a prescribed phase, using LS criterion for approximation. Several examples show the efficiency of the presented algorithm and its application in twodimensional cases. Except a difficulty encountered with this technique is high computation complexity when the size of image increases.
2
In practical situation, the technique to be applied for reconstructing an (N1 xN2) point image from the phase spectrum of some unknown observed image is as follows:
- First take some records including (L1 xL2) > 2(N1 xN2) samples of the phase spectrum at radians frequencies equally spaced between 0 and n .
- Arranging the phase samples as 1-D sequence.
- Apply the algorithm to the reconstructed sequence of each value of (L x L2) .
- Rearranging the 1-D reconstructed sequence as 2-D sequence (image).
- Use a direct comparison by using the 2-D correlation between the original and reconstructed image, this will identify the quality of reconstruction.
REFERENCES
[01] T. G. Stockham, T. M. Cannon, and R. B. Ingebreston. Blind deconvolution through digital signal processing. Proc. IEEE, pp. 678-692, Apr. 1975.
[02] H. C. Andrews and B. R. Hunt, Digital Image Restoration. Englewood Cliffs, NJ: Prentice Hall, 1977.
[03] S. R. Curtis, A. V. Oppenheim, J. S. Lim. Signal reconstruction from Fourier transform sign information, IEEE Trans. Acoust. Speech Signal Process. Assp-33 (1985), pp. 643-657.
[04] M. H. Hayes, J. S. Lim, A. V. Oppenheim. Signal reconstruction from phase or magnitude. IEEE Trans. Acoust. Speech Signal Process. Assp-28 (1980), pp. 672-680.
[05] G. A. Merchant, T. W. Parks. Reconstruction of signals from phase: Efficient algorithms, segmentation, and generalizations. IEEE Trans. Acoust. Speech Signal Process. Assp-31 (1983), pp. 1135-1147.
[06] H. Pozidis, A. P. Petropolu. Signal reconstruction from phase only information and application to blind system estimation. IEEE Internat. Conf. Acoust. Speech Signal Process. (ICASSP’97), Munich, Germany, 1997.
[07] T. F. Quatieri, A. V. Oppenheim. Iterative techniques for minimum phase signal reconstruction from phase or magnitude. IEEE Trans. Acoust. Speech Signal Process. Assp-29 (1981), pp. 1187-1193.
[08] B. Yegnanarayana, T. R. Krishman. Significance of group delay functions in signal reconstruction from spectral magnitude or phase, IEEE Trans. Acoust. Speech Signal Process. Assp-32 (1984), pp. 610-623.
[09] B. Yegnanarayana, S. T. Fatima, H. Murthy. Reconstruction from Fourier transform phase with applications to speech analysis. in IEEE Internat. Conf. Acoustic. Speech Signal Process. (ICASSP’87), Dallas, 1987.
[10] N. S. Hamid, T. F. Quatieri, J. S. Lim. Signal reconstruction from short Fourier time Fourier transform magnitude. IEEE Trans. Acoust. Speech Signal Process. Assp-31 (1983), pp. 986-998.
[11] A. P. Petropolu, H. Pozidis. Phase reconstruction from bispectrum slices. IEEE Trans. Signal Process. 46 (1998), pp. 527-531.
[12] H. Pozidis, C. L. Nikias. Signal reconstruction from the phase of the bispectrum. IEEE Trans. Signal Process. 40 (1992), pp. 601-610.
[13] T. Q. Nguyen, T. L. Laakso, D. Koilpillali. Eigenfilter approach for the design of all pass filters approximating a given phase response. IEEE Trans. Signal Process. 42 (1994), pp.2257-2257.
[14] J. M. Rouvaen, A. Djebbari, A. Menhaj-Rivenq, Al. Djebbari and Z. Chama. An algorithm for signal reconstruction from its phase based upon least square criteria. Signal Processing 79 (1999), pp. 67-72.
[15] J. Franklin. Matrix Theory. Prentice Hall, Englewood Cliffs, NJ, 1968.
