Recursive two-dimensional spectral factorization

CmCUITSSYSTEMSSIGNAlPROCESS VOL, 3, [NO, 2, I9~4

RECURSIVE Two-DIMENSIONAL SPECTRAL FACTOR!ZATION* L u & F. Chaparro ~

Abstract. The factorization of two-dimensional spectral functions is a problem of interest in the digital processing of signals because of the nonexistence of a fundamental theorem of algebra for bivariate polynomials. In this paper we consider the special case of rational spectral functions, which are quotients of polynomials or pseudopolynomials of finite order. We first characterize the general use of asymmetric halfplane factors, and then propose recursive procedures to obtain the factors. In the case of strictly positive spectral functions, we use the homomorphic transform to characterize the causal factor by its first-order data and then apply the modified leastsquares approximation to obtain the rational factors. For those cases of partially characterized spectral functions, we propose extensions of one-dimensional procedures. These alternative procedures are based on the modified least-squares approximation as well as on the two-dimensional Levinson algorithm. All the procedures presented are recursive, thus computationally efficient, and guarantee the stability of the factors in most cases. 1. Introduction Two-dimensional spectral factorization is a problem o f great interest in the area o f digital signal processing. Such factorization is needed, for instance, in the restoration o f images by Wiener [1] and Kalman [2] filtering and in the design and stabilization o f two-dimensional recursive filters [3,4]. The difficulty in performing such factorization lies in the fact that there is no fundamental theorem o f algebra for bivariate polynomials. O u r objective in this paper is to investigate the class o f two-dimensional rational spectral functions for which a rational factorization exists. These spectral factors are quotients o f polynomials or pseudopolynomials o f finite order and must be stable. A l t h o u g h it is possible to obtain these spectral factors by h o m o m o r p h i c transformation [5], they are usually o f infinite length. * Received September 14, 1983; revised December 7, 1983. This work was supported by NSF Grant ECS-8105988. ' Department of Electrical Engineering, University of Pittsburgh, Pittsburgh, Pennsylvania 15261, USA.

194

CHAPARR0

To obtain stable and finite-length factors, complicated nonlinear procedures are needed [6]. In this paper we extend the results of Murray [7], for the quarter-plane factors case, to the more general case of asymmetric half-plane factors. These results will clarify that, as indicated by Bose [8] and Woods [9], there are some rational spectral functions which cannot be realized as the amplitude response of a stable recursive filter. We then propose recursive procedures to obtain the rational spectral factors. Using the two-dimensional homomorphic transformation, we characterize the causal factor of the spectral function by its impulse response, and obtain from this data the rational factors using the modified least-squares (MLS) approximation procedures [10]. The parameters of the factors are obtained by solving a set of linear equations, which can be done recursively by relating this problem to the theory of orthogonal polynomials on the unit bidisk. For those cases for which the cepstrum cannot be calculated or for which there exists only a partial characterization of the spectral function, we propose some alternative procedures. These procedures are extensions to two-dimensions of those of Rissanen and Kailath [11], Bauer [12], and Scharf and Luby [13]. These procedures use the MLS approximation and the two-dimensional Levinson [14] algorithms. We will especially discuss the connection that exists among the two-dimensional Bauer algorithm and the theory of orthogonal polynomials and the Cholesky matrix factorization. To illustrate the various procedures as well as the issues involved, we provide some computer examples at the end of this paper. 2. Main results

A spectral function S (zl ,z2 ) is the z-transform of an autocorrelation or nonnegative-definite bisequence {r ( m, n ) 1. In this paper, we consider only spectral functions which are rational or have been approximated rationally. By rational spectral factorization we mean the decomposition of S (zl,z2) into factors which are the quotients of finite-order polynomials or pseudopolynomials, and that have no poles or zeros in a certain region of the space of two-dimensional complex variables. We call the polynomial case the quarter-plane (QP) factorization and the pseudopolynomial case the asymmetric half-plane (AHP) factorization. A. CEPSTRALPROCEDURE. According to [7], if the given rational spectral function S (zl,z2) is positive on the distinguished boundary, and if it admits a quarter-plane factorization, then the Fourier coefficients of log S (zl,z2), /z,I = 1, i = 1,2, are zero for all pairs of integers (m,n) in the second and fourth quadrants of its support plane. Also, the causal factor has no poles or zeros outside of the unit bidisk. For the more general case of asymmetric half-plane factors, we can obtain an analogous result in the following manner. Suppose that S (zl, z2 ) is positive

RECURSIVE Two-DIMENSIONAL SPECTRAL FACTORIZATION

195

on the distinguished b o u n d a r y and is factorized as S(z,,z2)

= H ( z ~ , z 2 )'

H"

tz~-1 ,z2-1 )

(1)

where H(z~,z2) -

(2)

B(z~,z2)

A (z~,z2) and B (z~ ,z~), A ( z l , z 2 ) are pseudopolynomials with a support R = {(re, n)[ m > O,n + m N >_ 0} . Also, let us assume that H ( z l , z 2 )

(3)

has no poles or zeros in the d o m a i n

D = DI UD2 ,

(4a)

where DI = {(zl,z2)] ]z2] = 1, !z~l ~ 11 ,

(4b)

D2 : [(zJ,z2)[[z2[-> 1,[zl[ = o~} .

(4c)

Theorem. If the above spectral function, S ( z l ,z2 ), admits an A H P factorization, then the Fourier coefficients o f log S ( z l ,z2 ) are zero for pairs o f integers (re, n ) such that m ( m N + n ) < O. N

If we let z~ = w z 2 , then S ( z l , z 2 ) is transformed into a quarter-plane factorable flanction [15], that is,

Proof.

(5)

S ( w , z 2 ) = H ( w , z2) H(w-',z-~ ~) ,

where H ( w, z2 ) is a quarter-plane filter with support in the first q u a d r a n t o f the support plane and has no poles or zeros outside the unit bidisk. Then according to M u r r a y ' s result, log S ( w , zz) = ~

g ( k , l ) w - k z2- l ,

Iwl = Izz] = I ,

(6)

(k,l) -N

with s = 0 for k l < O. Now, if we first let w = z l z 2 and then use the following one-to-one and onto coordinate transformation,

Iml E101E t n

-N

1

e

(7)

196

CHAPARRO

we get logS(z,,z2) = G

(8)

Y ( m , n + m N ) z-~m z 7

(re,n)

with d(m,n + mN) = 0

for m ( n + N m ) < 0 .

Remarks 1. The converse of this theorem holds if S ( z ~ , z 2 ) is a product of finiteorder pseudopolynomials. 2. The support of log S (zl,z2) depends on the support of the A H P factors, but it is always the union of two sectors each with a vertex at ( 0 , 0 ) and with a n aperture angle of less than 7r. We now propose a recursive procedure to obtain the rational spectral factors. If S ( z l , z 2 ) can be factored as in Equation (1), the region of support of the A H P filter H ( z ~ , z 2 ) can be obtained from the region of support of log S ( z , , z 2 ) . Using this information we can then transform this filter into a quarter-plane filter by a one-to-one and onto transformation [15]. It is then clear that it is only necessary to consider the quarter-plane factorization and apply the inverse transform. Suppose that a rational spectral function S (z~ ,& ) has a QP factorization; then -1

-1

logS(z,,z=) = logH(z,,z2) + log H ( z , ,z2 ) ,

(9)

where H ( z , , z 2 ) is a minimum-phase QP filter. Therefore h ( m , n ) and /~ ( m, n ), the Fourier coefficients of H (Zl,Z2) and log H (z,, z2 ), respectively, have the same support and one can derive recursive formulas [16] to obtain h ( m, n ) from /z ( m, n ), and conversely. These formulas give h(O,O) = exp (/~(0,0)) ,

7

I2 k=O

h(m,n)

f h ( m - k , n - f) /~(k,f),

n~O,

f=O

(10)

= t

I

k h ( m - k , n - f) /~(k,f), k=O

m~O.

,~

Now, given the characterization of H (z, ,z2) by its impulse response and

197

RECURSIVE Two-DIMENSIONAL SPECTRAL FACTORIZATION

the information that this filter is a stable QP filter, the problem is how to obtain the parameters of the approximating rational filter that matches the data. We propose to solve this problem by the modified least-squares approximation procedure [10], which consists of finding parameters [a ( m, n ), b ( m, n ) ] of a filter: H(Zl,Z~)

-

B(z~,z~)

(11)

~(z~,z~)

that minimize the error criterion, 1

('2 re~2 re J J IH(z~,z~) f t ( z , , z ~ ) -

e=-

471"2 0

B(z,,zz)pdw, dw2,

0

i = 1,2 .

z, = e I~',

(12)

This error measure is a modification of the canonical least-squares error measure. The advantages of using c in our case are: the resulting normal equations are a set of linear equations because e is quadratic with respect to the parameters of H(z~,z~), and only a finite characterization of H ( z , , z 2 ) by its impulse response and autocorrelation values is needed to get H ( z ~ , z ~ ) . Using the canonical least-squares error, the minimization is a nonlinear programming problem, which is computationally expensive, and H ( z l , z 2 ) needs to be characterized by the entire impulse response. From Equation (12) we obtain the following normal equations: (s,t) eSo ,

a(k,g)Ks.k(t,O = c~(M,N) 5(s,O)b(t,O),

(13a)

(k,e)~s

b(s,t)

a(k,Oh(s-k,t-O,

(s,t) E S~ ,

(13b)

(k,e)~So

where S~, Sb are the support lattices of A (z~ ,z2 ), B (zl ,z2 ) respectively, K,.k(t,e) = r ( s - k , t - O

-

~

h(i-k,j-e)h(i-s,j-t),

(i,j)~S b

(k,f), ( s , t ) E S o ,

(14)

and oe ( M , N ) is the non-negative minimum error for an A ( z l , z 2 ) of order ( M , N ) . Polynomials are ordered according to their highest lexicographic order monomial. Without loss of generality, we let A (z~,z2) and B (zl ,z2 ) have the same order, and the matrix form of the normal equations is then b(M,N) = H(M,N)a(M,N)

,

K ( M , N ) a ( M , N ) = ce(M,N) [1 0 ... 0V ,

(15a) (15b)

198

CHAPARRO

where b ( M , N ) , a ( M , N ) are vectors with entries the parameters {b ( k , f ) , The entries of the symmetric matrix K ( M , N ) are those defined in Equation (14), and this matrix can be expressed as

a (k,f)}.

= P(M,N) - H r ( M , N )

K(M,N)

where P ( M , N ) ,

=

P0

F~

r(i,-l)

(16)

Pf

ro

. . . FM

(17a)

P~,

r(i,O)

,

are the block Toeplitz matrices

H(M,N)

r (M,N)

r

H(M,N)

Po

r(i, 1)

r(i,N)

r(i,O)

Y,=

O M, I > N, and the matrices t~k,e

i - j ) ] , 0 _< i , j 0

for Izll = 1, i = 1,2. W e t o o k L = I = 7, calculated the matrix BL in E q u a tion (38), and the last row o f this matrix gave us the factor C(z~,z2)

= 1.002 - 0.3720 z~~ - 0.4703 z~~ + 0.2130 z~~ z~~

The o~ ( m , n ) decreased monotonically in both directions and remained unchanged for the high-order polynomials, indicating convergence o f the algorithm. W e then let V(z~,z2)

= I(1 - z ~ ) ( 1 - z~l)(1 + z-~tz-21)l 2 >- O,

Izll

= 1,

i = 1,2 ,

RECURSIVE Two-DIMENSIONALSPECTRAL FACTORIZATION 205

and took L = I = 7. In this case the c~( m , n ) behavior was completely different than in the previous case; the ~'s did not decrease monotonically and had similar values for low- as well as for high-order polynomials. We took this as a sign of no convergence. 4. Conclusions In this paper we have explored the characterization of two-dimensional spectral functions for the problem of spectral factorization, and have proposed recursive procedures to obtain the factors. These procedures are computationally efficient because of their recursive nature and simple implementation. We remark that there does not seem to be any difficulty in extending the results in this paper to higher dimensions. With respect to the stability of the spectral factors, we have found that in the cepstral procedure the minimum phase of the factors is directly related to the accuracy of the rational approximation. Since, at least theoretically, the modified least-squares approximation permits us to obtain a perfect match of the impulse response of the rational causal factor, it is therefore expected that the resulting approximating factors are minumum phase. In the alternate procedures, it is not possible to guarantee the stability of the spectral factors, but as shown iLn [10], stability occurs quite frequently. The two-dimensional Bauer algorithm behaves very much as in the onedimensional case. For the case when the spectral function is strictly positive on the distinguished boundary, the algorithm converges to a solution that is unique up to a factor of unit magnitude. Also, the connection between the existence of this solution and the necessary condition for factorization becomes obvious. In this case, when we obtain the factors using orthogonal polynomials, we find that the normalizing values decrease monotonically as the order of the polynomial increases. On the other hand, when the spectral function is zero for some fl'equencies on the distinguished boundary, the normalizing values do not decrease monotonically as the order of the polynomial increases, and they might even become zero or negative. Finally, it is of interest to note that the proposed procedures in this paper can be used as stabilization procedures for two-dimensional recursive filters, for the cases when the denominator of the unstable filter has no singularities on the distinguished boundary. The necessary conditions for rational factorization determine the cases for which the stabilization is not possible. Either the cepstral or the Bauer factorization can be used to factorize the magnitude square denominator into minimum-phase factors. These procedures are, for this particular case, more efficient than the Shanks stabilization procedure, since the factorization is done directly instead of in two steps as in the Shanks procedure, and also because of the ability to determine a p r i o r i the feasibility of the factorization which cannot be done in the Shanks procedure.

Acknowledgment The author wishes to thank the reviewers for their constructive criticism.

206

C~APARRO

References [1]

M . P . Ekstrom, "Realizable Wiener filtering in two dimensions," IEEE Trans. ASSP, Vot 30, February 1982. [2] J.W. Woods, "Markov image modeling," IEEE Trans. AC, Vol. 23, October 1978. [3] J.L. Shanks et al., Stability and synthesis of two-dimensional digital filters," IEEE Trans. AEA, Vol. 20, June 1972. [41 Y.V. Genin and Y. G. Kamp, "Two-dimensional stability and orthogonal polynomials on the hypercircle," Proc. IEEE, Vol. 64, June 1977. [5] M . P . Ekstrom and J. W. Woods, "Two-dimensional spectral factorization with applications in recursive digital filtering," IEEE Trans. ASSP, Vol. 24, April 1976. [6] M . P . Ekstrom et al., "Two-dimensional recursive filter design - a spectral factorization approach," IEEE Trans. ASSP, Vol. 28, February 1980. [71 J.J. Murray, "Spectral factorizafion and quarter-plane digital filters," IEEE Trans. CAS, Vol. 25, August t978. [8] N. K. Bose, "Problems on stabilization of multidimensional filters via Hilbert transforms," IEEE Trans. GE Vol. 12, October 1974. [9] J . W . Woods, IEEE Trans. GE (Corresp.), Vol. 12, July 1974. [10] L. F. Chaparro and E. I. Jury, "Rational approximation of 2-D linear discrete systems," IEEE Trans. ASSP, Vol. 30, October 1982. [11] J. Rissanen and T. Kailath, "Partial realization of random systems," Automatica, Vol. 8, 1972. [12] D. C. Youla and N. N. Kazanjian, "Bauer-type factorization of positive matrices and the theory of matrix polynomials orthogonal on the unit circle," IEEE Trans. CAS, Vol. 25, February 1978. [13] L. L. Scharf and J. C. Luby, "Statistical design of autoregressive moving average digital filters," IEEE Trans. ASSP, Vol. 27, June 1979. [141 J. H. Justice, "A Levinson-type algorithm for two-dimensional Wiener filtering using bivariate Szeg6 polynomials," Proc. IEEE. Vol. 65, June 1977. [151 B.T. O'Connor and T. S. Huang. "Stability of two-dimensional recursive filters," IEEE Trans. ASSP, Vol. 26, December 1978. [16] D. E. Dudgeon, "The computation of two-dimensional cepstra," IEEE Trans. ASSP, Vol. 25, December 1977. [171 F. R. Gantmacher, The Theory of Matrices, Vol. 1, Chelsea, New York, 1960.