IEEE TRANSACTIONS O N ACOUSTICS, SPEECH, SIGNAL AND

PROCESSING, VOL. ASSP-29, NO. 6 , DECEMBER 1981

1187

Iterative Techniques for Minimum Phase Signal Reconstruction from Phase or Magnitude

Absstract-In thispaper, wedevelop iterativealgorithmsforreconstructing a minimum phase sequence from phase or magnitude ofits the Fourier transform. These iterative solutions involve repeatedly imposing a causality constraint in the time domain and incorporating the known phase or magnitude function in the frequency domain. This approach i the basis of a new means of computing the Hilbert transform s of the logmagnitude or phase of the Fourier transform of a minimum phase sequence which does not require phase unwrapping. Finally, we discuss the potential use of this iterative computation in determining samples of the unwrapped phaseof a mixed phase sequence.

quency domain. Another iteration in this same style recovers a finite length mixed phase signal from the phase of its Fourier transform by imposing a finite length constraint in the time domain and the known phase in the frequency domain [4]. In this paper, we begin in Section I1 with a discussion of a number of equivalent conditions for a sequence to be minimum phase. In Sections I11 and lV,we use these conditions in developing two iterative reconstruction algorithms for minimum phase signals, one for reconstruction when the phase is known and the other for reconstruction when the magnitude is I. INTRODUCTION known. NDER certain conditions a signal can be reconstructed In Section V, we discuss the discrete Fourier transform from a partial specification in the time domain, in the (DFT) realizations of the algorithms and illustrate the reconfrequency domain, or in both domains. A minimum or maxi- struction process with examples. mum phase signal, in particular, can be recovered from the In Section VI, we propose the use of the algorithms of Secphase or magnitude of its Fourier transform [l] . The conven- tions I11 and I in implementing the Hilbert transform. Of V tional reconstruction algorithm involves applying the Hilbert particular importance is reconstruction of the log-magnitude transform to the log-magnitude or phase of the Fourier trans- from phase since the proposed iterative approach requires only form to obtain the unknown component. the principal value of the phase, while the direct DFT impleIn this paper, we take an alternative approach by developing mentation of the Hilbert transform requires the unwrapped iterative algorithms for reconstructing a minimum (or maxi- phase [5]. The proposed technique,therefore, avoids probmum) phase signal from the phase or magnitude of its Fourier lems typical of phase unwrapping such as detection of the distransform. Specifically, we develop algorithms which impose continuities in the principal value of the phase [l], [6] . Also, causality in the time domain and the given phase or magnitude in Section VI, we suggest the use of this new approach to imin the frequency domain, in an iterative fashion. plementing the Hilbert transform as the basis for a phase unIterative algorithms similar to those we discuss here have wrapping algorithm. been useful in a number of areas where partial information in thetwo domains is available. Inparticular,the algorithms 11. THE MINIMUMPHASECONDITION presented in this paper are similar in style to the GerchbergIn general, a signal cannot be uniquely specified by only the Saxton algorithm [2] and an iterative algorithm by Fienup phase or magnitude of its Fourier transform. However, one [3] , in alternately incorporating partial information in the condition under which the magnitude and phase are related is time and frequency domains. The Gerchberg-Saxton algo- the minimum phase condition and under this condition a signal rithm recovers a two-dimensional complex signal by iteratively can be uniquely recovered from the magnitude of its Fourier imposing the finite extent of the signal in the space domain transform or to within a scale factor, from the phase of its and its magnitude in both the space and frequency domains. Fourier transform. In this section, we discuss anumber of Similarly, Fienup's algorithm recovers a real two-dimensional equivalent conditions for a signal to be minimum phase. These signal by iteratively imposing the finite extent and positivity conditions will be of particular importance in Section 111 in deof the signal in the space domain and its magnitude in the fre- veloping the iterative algorithms. In the following discussion we restrict the z transform of Manuscript received June 13,1980; revised June 30, 1981. Thiswork the sequence h(n) to be a rational function, which we express was supported in part by the Department of the Air Force under Conin the form tractF19628-80-C-0002,inpart by theNationalAeronauticsand
Space Administration under Grant NAGS-2, and in part by the Advanced Research Projects Agency, monitored the Office of Naval Reby search under Contract NOOO14-75-C-0951-NR 049-328. The views and conclusions contained in this document are those the contractor and of should not be interpreted as necessarily representing the official policies, either expressedor implied, of the US.Government. T. F. Quatieri, Jr. is with Lincoln Laboratory,Massachusetts Institute of Technology, Lexington, MA 02173. A. V. Oppenheim is with the Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA 02139.

M i

H(Z) = A Z n o

k=1

n n
=

MO

(1 - ~

k=1

P i

(1 - C k z - l )

k=l ldkl

rI (1 n (1
)

- bkz)

PO

(1)
- dkZ)

where la,/, Ibkl, Ickl, and

areless

thanor equal to

0096-3518/81/1200-1187$00.75 0 1981 IEEE

1188

IEEE TRANSACTIONS ACOUSTICS, ON SPEECH,

AND SIGNAL PROCESSING, VOL.

ASSP-29, NO. 6, DECEMBER 1981

unity, zno is a linear phase factor, and A is a scale factor. When, in addition, h(n) is stable, i.e., En Ih(n)l < 00, Ick and Idk\ are strictly less than one. A complex function H(z) of a complex variable z is defined to be minimum phase if it and its reciprocal H- (z) are both analytic for I > 1. A minimum phase sequence is then dez fined as a sequence whose z transform is minimum phase. For H(z) rational, as in (2), the minimum phase condition excludes poles or zeros on or outside the unit circle in the z plane or at infinity. As a consequence, the factors of the form (1 - bkz) corresponding to zeros on or outside the unit circle and the factors of the form (1 - dkz) corresponding to poles on or outside the unit circle will not be present. Furthermore, in (l), no = 0 to exclude poles or zeros at infinity. Thus, for H(z) minimum phase, (1) reduces to

Again, from (2) it follows that these conditions are necessary since (2) has no poles or zeros outside the unit circle or at infinity, guaranteeing causality, and from the initial value theorem h(0) = lim H(z) = A . To demonstrate that these conz+m ditions are sufficient, we note that again causality of h(n) will eliminate factors of the form (1 - d k z ) in the denominator of (1). Furthermore, since the conditions require that h(n) be causal, the initial value theorem can be applied with the result that

h(0) = lim H(z) = lim AZno

n
k=l

(1 -

Since h(0) = A ,
M,

M i

n (1
Pi

akz-)

H(z) = A

k=l

k=l

rJ

(2)
(1 - c k z - l )

where lak] and ck are both strictly less than unity. From (2) other conditions can be formulated for a signal to be minimum phase. Two conditionsin particular which we discuss below are particularly useful in the context of the iterative algorithms to be discussed in Sections I11 and IV. Minimum Phase Condition A Consider h(n) stable and H(z) rational in the form of(1) with no zeros on the unit circle. A necessary and sufficient condition for h(n) to be minimum phase is that h(n) be causal, i.e., h(n)= 0, n < 0, and no in (1) be zero. From (2), it follows that these conditions are necessary. To show that they are sufficient, we want to show that they force (1) to reduce to (2). Clearly, factors of the form (1 - d k z ) , ldkl < 1 in the denominator introduce poles outside the unit circle which would violate the causality condition since h(n) is restricted to be stable. With no = 0 in (l), factors of the form (1 - bkz) would introduce positive powers of z in the Laurent expansion of H(z), requiring h(n) to have some nonzero values for negative values of n , thereby again violating the causality condition. Therefore, these factors cannot be present and with no = 0, (1) reduces to (2). Finally, because our condition assumes h(n) is stable and that H(z) has no zeros on the unit circle, h(n)is minimum phase. The above minimum phase conditions require that h(n) be causal and that the unwrapped phase function have no linear phase component. Another slightly different set of necessary and sufficient conditions for a signal to be minimum phase can be stated as follows. Minimum Phase Condition B Consider h(n) stable and H(z) rational i the form of (1) n with no zeros on the unit circle. A necessary and sufficient condition for h(n) to be minimum phase is that h(n) be causal, Le., h(n) = 0, n < 0 and h(0) = A where A is the scale factor of (1 1.

and since Ibk < 1 this requires that no = 0 and the bks be equal to zero. Thus, again (1) reduces to (2). Another condition which can be shown to be equivalent to minimum phase condition A or B or our original definition of a minimum phase sequence is that the log-magnitude and unwrapped phase of H ( a ) are related through the Hilbert transform [ l ] . The Hilbert transform relation guarantees thata minimum phase sequence can be uniquely specified from the Fourier transform magnitude and, to within a scale factor, from the Fourier transform phase. One technique for minimum phase signal reconstruction from phase or magnitude relies on a DFT implementation of the Hilbert transform [5]. In the next two sections, we take an alternate approach which invokes an iterative computation. Motivated by the minimum phase condition A, when the phase is given weimpose, in an iterative fashion, causality in the time domain and the known phase in the frequency domain. When the resulting sequence satisfies minimum phase condition A and has the given phase, it must equal h(n) to within a scale factor. Likewise, motivated by the minimum phase condition B, when the magnitude isgiven,we impose, in an iterative fashion, causality and the initial value h(0) in the time domain, and the known magnitude in the frequency domain. When the algorithm results in a sequence which satisfies minimum phase condition B and has the given magnitude, it must equal h(n).
111. AN ITERATIVE ALGORITHM SIGNAL FOR RECONSTRUCTION FROM PHASE The iterative algorithm for reconstructing a minimum phase signal from its phase function is shown in Fig. 1. The function hk(n) represents the signal estimate on the kth iteration and &+ (n) = h&) u(n) where u(n) is theunitstepfunction. The function e,(u) is the known phase and M k + ( w ) and r3k+l (a)are the Fourier transform magnitude and phase of h;c+ (n),respectively. The algorithm begins with an initial guess Mo(w) of the desired Fourier transform magnitude and the inverse Fourier transform of Mo(u) [jeh(u)] is taken. This step yields exp h,(n), the initial estimate of h(n). Next, causality is imposed so that ho(n) is set to zero for n < 0 to obtain g,(n). The

1.

QUATIERI PHASE AND OPPENHEIM: MINIMUM

SIGNAL RECONSTRUCTION

1189

Therefore, from (9) and (lo), and the identity lexp [jeh(o)] = 1:

1'

- Mk + 1(a) exP [jeh(o)112 d o

Fig. 1. Iterative algorithm to recover h(n) from its phase.

2g
phase of the Fourier transform of i l ( n ) is then replaced by the given phase and the procedure is repeated. We now show that the mean-square error between h(n) and hk(n) or, equivalently, between their respective Fourier transforms, H ( o ) = Mh(o) [je, ( ) and &(o) =Mk(w) exp exp a] [jOh(w)] ,is nonincreasing on successive iterations. The meansquare error on the kth iteration from Parseval's Theorem can be written as

1,
n

IH(o)-Hk+l(o)1' d o

> E k + 1*

(1 1)

Since Ek is, therefore, nonincreasing and has a lower bound of zero, Ek must converge to a unique limit [7] . The nonincreasing nature of E,, however, is not sufficient to guarantee that the iterates hk(n) converge. Nevertheless, V a converging solution with a rational z transform exists, we can show that

where a! is a positive constant. To see this, note from (6), (7), and (10) that the equality in (1 1) holds if and Only if hk(n) = h;, + 1(n) = 0 for n < 0, and e h ( o )= ek + (0). Therefore, since O,(w) contains no linear phase component (i.e., no = 0), if hk(n) converges to a sequence whose z transform is of the form in (l), the converging solution must satisfy the minimum phase condition A. Consequently, the converging solution is minimum phase with phase e h ( o ) , and (12) must hold.' When h(n) is of finite duration (i.e., H(z) has no poles), we can impose not only causality, but also a finite duration constraint within the iteration. Under these particular constraints, the DFT realization of our iterative procedure (see Section V) always converges to a limit of the form in (12) [8].

"x
n

Ih(n)- ik+l(n)['*

Iv. AN ITERATIVE ALGORITHM SIGNAL FOR RECONSTRUCTION FROM MAGNITUDE

Next, from Parseval's Theorem, we write (8) as

&

-Mk+i(a)exp

[ i e k + i ( ~ ) l 1 2d o .

With the triangle inequality for vector differences, we have at

'The constant 01 in (12) is constrained to be positive since a negative value introduces an additive factor of II into the phase function.

1190

IEEE TRANSACTIONS ON ACOUSTICS, SPEECH, AND

SIGNAL PROCESSING,VOL.ASSP-29,NO.

6 , DECEMBER 1981

F- I

hk(nl

IMPOSE CAUSALITY
AND h (01

-t + l ( n ) = h
I
F

hk(n)

n > 0

h(O1 0

n = 0
n c O

~ ( w+) e x p [,ek+, ( W I ] ~

Mk+l(~I--+Mh(~)

~,(w~exp[je,+,(w)]

Fig. 2. Iterative algorithm to recover h(n) from its magnitude.

The function M h ( a ) is the known magnitude and M k + (a) and Ok + (a) the Fourier transform magnitude and phase are " of & + (n), respectively. > 27l -" p " ) - & + 1 ( 4 1 2 ' dm The algorithm begins with an initial guess Oo(a) the deof sired phase, and the inverse transform of Mh(a) [jeo(a)] exp 2Ek+1(2 1) is taken. This step yields ho(n), the initial estimate of h(n). Next, on the basis of the minimum phase condition B, causalSince Ek is nonincreasing and has a lower bound of zero, it ity and the known value of h(0) are imposed so that ho(n) is must converge to a limit point [7]. set to zero for n < 0 and set to h(0) for n = 0, to obtain il (n). As withthe algorithm in Section 111, although we have The magnitude of the Fourier transform of i l ( n ) is then reshown that the error Ek is nonincreasing, we have not shown placed by the given magnitude and the procedure is repeated. that the iterates hk(n) converge. However, ifthe iterates conIt has not been possible to show that the mean-square error, verge to a sequence whose z transform is rational with no zeros as considered in Section 111, is nonincreasing forthis algoon the unit circle and which is causal with initial value h(O), rithm. However, an error function that is nonincreasing is the from the minimum phase condition B, the converging solution mean-square difference between the known magnitude and the must be minimum phase. Consequently, if in additionthe estimate M k ( a ) on each iteration, i.e., magnitude of the Fourier transform of the converging solution equals Mh(w), the solution is the unique minimum phase sequence associated withMh(a), i.e., h(n). The convergence of hk(n) has yet to be rigorously proven even when a finite length constraint is imposed within the iterTo show that Ek is nonincreasing, we first use the identity ation [8]. Empirically, however, we have found the DFT reallexp [ j O k ( a ) ] = 1 to express Ek as ization of the algorithm to always converge. In the next section, we shall illustrate the convergence of h,(n) to h(n) with 1 " an example. Ek = I'h(a> - M k ( a ) 1 2 exp [iek(a)lI2

1'

4,

From Parseval's Theorem, (15) is givenin the time domain by

Ek =
n

Ihk(n) - $k(n)l2.

(16)

From (1 3), it follows that

OF ALGORITHMS V. REALIZATIONS THE ITERATIVE USING THE DFT Since the iterative algorithms will be implemented on a digital computer, we can compute a Fourier transform at only a finite number of points. In particular, we shall use the DFT. One consequence of the DFT realization is that our desired sequence h(n) must be of finite duration. Imposing a finite duration constraint within the iterations, however, does not change the nonincreasing nature of the error functions, as can be seen from (8) and (19). A second consequence of the DFT realization is that only uniformly spaced samples of the phase and magnitude functions are available. Nevertheless, it is again possible to show

QUATIERI AND OPPENHEIM: MINIMUM PHASE SIGNAL RECONSTRUCTION

1191

128

256
4

384

Fig. 3. Convergence of h&) in example 1: (a) original, (b) 1 iteration, (c) 5 iterations, (d) 45 iterations.

1000

2000 3000 ------cHZ

4000

that the nonincreasing natureof Ek is not altered when we use samples of thl: magnitudes and Fourier transforms in the expressions for E , in (5) and (14) .[9] , [ 101 . Finally, questionsof convergence need to be addressed. Consider first, minirhum phak reconstruction from phase samples. When H(z) is constrained to have no conjugate reciprocal Zero pairs and no zeros on the. unit circle, a unique sequence h(n) of length M (to within a scale factor) is guaranteed when given M - 1 or more phase samples of Bh(w) in the open frequency interval (0, T ) [9] . A minimum phase sequence, in particular, satisfies these constraints. Therefore, the DFT realization of the iterative algorithm to reconstruct a minimum phase sequence from its phasesamples can be implemented with a DFT of length N > 2M. Furthermore, this iteration will converge to ah(n) for 0 <n < N - 1 where a is positive 181 Consider next, the dual problem of developing a DFT realization of the iterative algorithm to recover a minimum phase sequence of lengthM from a magnitude function. In this case, there exists only one M point sequence, i.e., the minimum phase sequence h(n) when h(0) is specified along with M or more uniformly spaced samples of the magnitude in the halfopen frequency interval [0, [lo]. Therefore, a DFTrealization of the iteration can be implemented with DFT length N2 2M - 1. If the algorithm converges to a causal sequence of length M with initial value h(0) and the known magnitude samples, the converging solution must equal h(n) for 0 < n < N - 1. To iilustrate, we now consider two examples where DFT the length is 512 points, which is twice the length of h(n). In the first example, the initial magnitude guess is unity, and in the second example the initial phase guess is zero.

Fig. 4. Convergence of log IHk(w)l (in decibels) in example 1: (a) original, (b) 1 iteration, (c) 5 iterations, (d) 45 iterations.

I
0

i
O I P

0 7*&d L
c

.)

Fig. 5. Convergence of hk(n) in example 2: (a) original, (b) 1 iteration, (c) 5 iterations, (d) 25 iterations.

h(n). The functions h,(n) and log [hfk(id)] are depictedin Figs. 3 and 4 along with the originals for k equal to 1,5, and 45. The signal hk@) (to within a multiplicative constant) and the spectrum log [Mk(w)] (to within an additive constant) are indistinguishable from the originals after 45 iterations. Example 2: Signal Reconstruction from Magnitude In this example, we consider the sequence of example 1, but where theFouriertransformmagnitude isgiven. Thefunctions hk(n) and Bk(w) are depicted in Figs. 5 and 6 with the originals for k equal to 1, 5, and25.Thefunctions hk(n) and 13,(w) are indistinguishable fromthe originals after 25 iterations. VI. A BASIS FOR IMPLEMENTATION HILBERT OF THE TRANSFORM PHASEUNWRAPPING AND In this section, we proposetwocomputationalalgorithms based on the procedures of Sections I11 and IV: 1) an iterative

Example 1: Signal Reconstncction from Phase Consider a 256-point minimum phase signal h(n) illustrated in Fig. 3.2 The phase is known and we wish to reconstruct
2The z transform of this signal consists of two complex pole pairs and onecomplexzero pair all withintheunit circle. For n > 256, h(n) has decayed to effectively zero.

1192 TRANSACTIONS ACOUSTICS, IEEE ON SPEECH, SIGNAL AND PROCESSING,

VOL. ASSP-29, NO. 6 , DECEMBER 1981

+
Hz

Fig. 6 . Convergenceof e k ( w ) (in radians) inexample 2: (a) original, (b) 1 iteration, (c) 5 iterations, (d) 25 iterations.

approach to computing the Hilbert transform, and 2) the potential use of l ) as the basis of a phase unwrapping algorithm. For a minimum phase signal, the log-magnitude and phase of the Fourier transform are related through the Hilbert transform and the direct implementation of the Hilbert transform using the DFT has been extensively investigated [5] . One disadvantage of this implementation is that in computing the logmagnitude from the phase, samples of the unwrapped phase are required and are often difficult to compute. An alternative to the direct implementation of the Hilbert transform exploits the iterative algorithms of Sections I11 and IV. When the phase is given, through the use of the algorithm in Section 111, ah(n) is first obtained from the phase and, in particular, does not require samples of the unwrapped phase. From ah@) the log-magnitude of orH(w), representing the Hilbert transform of the phase to within an additive factor is then computed. Furthermore, with a fixed DFT length, by increasitig the number of iterations we can come arbitrarily close to samples of the log-magnitude. The direct approach, on the other hand, requires an increase in the DFT length for an increase in accuracy [l ] , [6] . A similar procedure can, of course, be applied through the use of the iterative algorithm in Section IV to implement the Hilbert transform of a given log-magnitude function. If h(0) is not known a priori [recall (13)] ,it can be obtained (at least in theory) from the magnitude, although in practice h(0) can be computed only approximately [l] . However, it was found empirically thattheiterates always converge to h(n) when only causality is imposed in the time domain (i.e., h(0) is assumed unknown) and the initial phase e,(w) is set to zero. This indirect approach to computing the Hilbert transform suggests a potential alternative to available phase unwrapping algorithms [1] , [6] . Let O(w) denote the desired unwrapped phase of the Fourier transform H(w) and O,(w) its value modulo 2n Furthermore, that the linear phase Omponent of H(w), i.e., no in (l), is known. The proposed phase unwrapping algorithm proceeds as follows. 1) R~~~~~ the linear phase component to obtain the prinCiPal value of the phase of H(w) exp [-inowl
e ^ p W
7

2) Apply the iterative algorit+ of Section I11 with a causality constraint and with phase t9,(w) to obtaina minimum phase sequence hmp(n). 3) Compute log (Hmp(o)I where H m p ( o ) is the Fourier transform of hmp(n). 4 Apply the Hilbert transform to log IHmp(w)l to obtain ) the unwrapped phase function B(w) - now. 5) Add the linear phase component to obtainthe desired unwrapped phase. Of particular interest is step 2 which yields the same minimum phase sequence h,,(n) that would be obtained by a Hilbert transform of the unwrapped phase, but bypasses the need for phase unwrapping. This algorithm has performed successfully on a number of simple mixed-phase :sequences (i.e., two and three poles and/or zeros) when the linear phase component was known exactly. Furthermore, it yielded the correct phase function when poles and zeros were placed close to the unit circle. There are several potential difficulties in the use of this algorithm.First,the minimum phase sequence hmp(n) derived from the iteration is of infinite extent regardless of whether the original sequence h(n) is of finite duration [lo] . Therefore, a possible problem with aliasing arises. The DFT length must be sufficiently large so thatthe minimum phase sequence h,,(n) decays effectively to zero. In particular, when hmp(n)= 0 for n > M , the DFT length, from our discussion in Section V, should be at least 2M. One possible procedure for removing the effects of aliasing invokes the principal value of the phase t9,(w) in a style similar to that in [6]. In particular, we might consider adding a step 6 of the following form.3 6 ) At each frequency,subtract 2n from the unwrapped phase estimate until the result is between -n and n. Then add this multiple of 2n so found to the principal value of the original phase function Bp(w). This will give an unwrapped phase free of the aliasing errors introduced in the iterative process. A second potential difficulty is the requirement that the linear phase factor of H(z) be known. Often, a priori knowledge of such a factor is difficult to obtain. One means of obtaining a linear phase estimate is to numerically integrate the phase derivative [I] . The sensitivity4 of our proposed phase unwrapping algorithm to deviations from the true linear phase in such an estimate is an area which needs to be explored if a practical algorithm is to evolve. VII. SUMMARY CONCLUSIONS AND In this paper, we have developed iterative algorithms for reconstructing a minimum phase sequence from either the phase or magnitude of its Fourier transform. When the phaseis known,the mean-square error between the desired Fourier transform and its estimate was shown to be nonincreasing on successive iterations. Likewise, when the magnitude is given,
3The authors acknowledge with thanks this suggestion by Dr. R. W. Schafer. 4When the residual phase linear is negative (so that no in znonegais tive), a causalconverging solution of the form in (1) derived from step 2 of our proposed algorithm will be mixed phase. In this case, the number of zeros outside the unit circle cannot be greater than In,,/. When no is positive, no causal converging solution of the form in (1) can exist.

denoted by

QUATIERI AND OPPENHEIM:MINIMUMPHASE

SIGNAL RECONSTRUCTION

1193

on successive iterationsthe mean-square error between the known magnitude and its estimate is nonincreasing. In addition, we noted that convergence of the iteration with known phase samples (i.e., the DFT realization) has been demonstrated, but convergence of the iteration with magnitude samples has been observed only empirically. Finally, we suggested two computational algorithms based on the iterative procedures: 1) a new means of implementing the Hilbert transform which avoids the need of an unwrapped phase, and 2) a new procedure for phase unwrapping. The iterative algorithms, as presented, rely on exact knowledge of the magnitude, phase, and the initial value of the desired signal. Sensitivity to the inexactness of these quantities, to quantization noise, and other forms of degradation is not understood and is an important area for future research. In practice, we have found that the iterative algorithms converge sometimes slowly (e.g., after several hundred iterations) and sometimes quickly (e.g., after a few iterations). Consequently, determining rates of convergence in terms of characteristics of the minimum phase signal and initial magnitude or phase estimates, and methods of speeding up convergence need to be explored. Another area being considered is the interchange of the signal reconstruction problems. In particular, we have found empirically thatwhen Mh(o)and Oh(o) are interchanged through j log H(o), slightly modified version of the iterative a algorithm of Section 111, requiring a phase function, wiu recover h(n) from its magnitude. Likewise, when the phase is known, h(n) is recovered by a procedure similar to the iteration in Section IV which requires a magnitude function. These results have led to some interesting theoretical speculations aboutthe duality of the reconstruction problems and their iterative solutions.

[9] M. Hayes, J. S . Lim, and A . V. Oppenheim, Signal reconstruction from phase and magnitude, IEEE Trans. Acoust,, Speech, Signal Processing, vol. ASSP-28, Dec.1980. [ l o ] T.F. Quatieri, Phase Estimationwithapplication to speech analysis-synthesis, Sc.D. dissertation, Dep. Elec. Eng. and Computer Science, M.I.T., Cambridge, MA, Nov. 1979.

REFERENCES
A. V. Oppenheim and R. W. Schafer, Digital Signal Processing. Englewood Cliffs, NJ: Prentice-Hall, 1975. R. W. Gerchberg and W. 0. Saxton, A practical algorithm for the determination of phase from image and diffraction planepictures, Optik, vol. 35, p. 237, 1972. J. R. Fienup, Reconstruction of an object from the modulus of its Fourier transform, Opt. Lett.,vol. 3, p, 27, July 1978. A. V:Oppenheim, M. H. Hayes, and J. S . Lim, Iterative procedures for signal reconstruction from phase, in Proc. 1980 Znt. Opt. Comp. Conf,vol. SPIE-231, Apr. 1980, pp. 121-129. B. Gold, A. V. Oppenheim, andC. M. Rader, Theory and implementation of thediscreteHilberttransform,in Proc. Symp. Comput. Processing Commun., vol. 19. New York: Polytechnic Press, 1970. J. M. Tribolet, A new phase unwrappingalgorithm, ZEEE Trans. Acoust., Speech, Signal Processing, vol.ASSP-25,Apr. 1977. M. Rosenlight, Introduction to Analysis. Glenview, IL: Scott, Foresman, 1970. V. T. Tom, T . F. Quatieri, M. H. Hayes, and J. H. McClellan, nonexpansive signal reconstruction Convergence of iterative algorithms, Lincoln Lab.,M.I.T., Lexington, MA, Tech. Note 1980-21, 1980. Also IEEE Trans. Acoust., Speech, Signal Processing, to be published.