Theory and Design of Signal–Adapted FIR Paraunitary Filter Banks∗

Theory and Design of Signal–Adapted FIR Paraunitary Filter Banks∗ Pierre Moulin†‡ and M. Kıvan¸c Mıh¸cak U. of Illinois, Beckman Inst. and ECE Dept, 405 N. Mathews Ave., Urbana, IL 61801 submitted to IEEE Trans. on Signal Processing, January 30, 1997. Revised, August 25, 1997. EDICS Number: SP 2.4

Abstract We study the design of signal–adapted FIR paraunitary filter banks, using energy compaction as the adaptation criterion. We present some important properties that globally optimal solutions to this optimization problem satisfy. In particular, we show that the optimal filters in the first channel of the filter bank are spectral factors of the solution to a linear semi–infinite programming (SIP) problem. The remaining filters are related to the first through a matrix eigenvector decomposition. We discuss uniqueness and sensitivity issues. The SIP problem is solved using a discretization method and a standard simplex agorithm. We also show how regularity constraints may be incorporated into the design problem so as to obtain globally optimal (in the energy compaction sense) filter banks with specified regularity. We also consider a problem in which the polyphase matrix implementation of the filter bank is constrained to be DCT–based. Such constraints may also be incorporated into our optimization algorithm, so we are able to obtain globally optimal filter banks subject to regularity and/or computational complexity constraints. Numerous experiments are presented to illustrate the main features that distinguish adapted and nonadapted filters, as well as the effects of the various constraints. The conjecture that energy compaction and coding gain optimization are equivalent design criteria, is shown not to hold for FIR filter banks.

∗

Work partly presented at ICASSP’95 and ICASSP’96.

†

Part of this work was performed while P. Moulin was at Bellcore, Morristown, NJ.

‡

Corresponding Author: tel (217) 244-8366, fax (217) 244-8371, Email: [email protected]

1

1

Introduction

One may broadly identify three main approaches to filter bank design. The first, classical approach consists in designing frequency–selective filters intended for applications such as speech coding and communications [1]. The second approach, motivated by developments in wavelet theory, consists in designing regular filters, and finds applications in subband systems using cascaded filter banks, e.g., in image coding [2]. The third approach, which forms the central topic of this paper, consists in designing signal–adapted filters, using energy compaction as the adaptation criterion [3, 4]. The use of energy compaction filters is motivated by applications in signal analysis, denoising, compression, and progressive transmission. The goal is to design the filter bank so that good L2 approximations to the signal of interest can be computed from a reduced number of channels. The case of 2–channel filter banks has received considerable attention, as illustrated by recent theoretical and experimental work in [5, 6, 7]. While only minor improvements over classical designs have been found for low–frequency signals, other types of of signal, such as image textures, greatly benefit from the adapted filter bank approach. In this paper, we present new theoretical developments and optimal design methods for M –channel FIR paraunitary filter banks (M ≥ 2). Such filter banks offer further advantages over adapted 2–channel filter banks, including a greater capability for isolating narrowband, high–frequency components of a signal [8]. As an elementary example, a narrowband (k+1)π signal component with 20% spectral occupancy (not necessarily in a frequency interval [ kπ ], 5 , 5

0 ≤ k < 5) would be conveniently extracted using an adapted 5–channel filter bank. Consider the M –channel FIR paraunitary filter bank in Fig. 1, which decomposes the input signal x into M components (M output channels). All filter coefficients are real–valued. While traditional (nonadapted) designs aim at optimizing some function of the filters, e.g., a weighted combination of stopband attenuations [1, Eq. (6.5.5)], the signal–adapted filter bank design problem considered in this paper is formulated as a multistage optimization problem [3, 4]: (EC) Design the components of the filter bank so that the energy in the first channel is maximized. Then use the remaining degrees of freedom to maximize the energy in the second channel, and repeat this operation successively for all remaining channels. It may not be obvious that there is any “degree of freedom” left after each step. However, a result due to Vaidyanathan et al. [9] shows that this is the case.

2

M. Unser has also conjectured that (EC) is equivalent to maximization of the coding gain 1 of the filter bank. However, the proof of this conjecture in [4] holds only for infinite–length filters (N = ∞) or for block transforms. In these two special cases, the solution respectively produces brickwall filters, and the Karhunen–Lo`eve transform (KLT) [3, 4]. For constrained–length filters, (EC) is a seemingly intractable nonlinear optimization problem for which no solution is known, and no algorithm has yet been found to be satisfactory [3, 4]. In this paper, we show how a globally optimal solution of (EC) may also be computed for constrained–length filters, extending a technique recently developed for the 2–channel case [7]. The polyphase representation of the filter bank plays a fundamental role in the solution to the optimization problem. We also show that (EC) and coding gain optimization produce different solutions in general. The remainder of this paper is organized as follows. In Sec. 2, we present the basic assumptions and the underlying optimization problems which arise from deterministic or stochastic formulations of the signal–adapted design problem. We also present an overview of the solution to the optimization problem. The difficult part of (EC) is the design of the first filter H0 (z), and we show how to solve that problem in Sec. 3. The key is the formulation of the latter optimization problem as a linear semi–infinite program (SIP) [7, 10] that may be solved using simple, globally optimal, and numerically stable algorithms. In Sec. 4, we present a characterization of the optimal filters in terms of their zeroes on the unit circle and discuss uniqueness and sensitivity issues. In Sec. 5, we describe a discretization algorithm for solving the SIP problem. In Sec. 6, we address the signal–adapted filter design problem under additional constraints, namely computational efficiency, and regularity constraints, which produce M -band wavelet filters. Numerical results are presented in Sec. 7, and our conclusions are given in Sec. 8.

1.1

Notation

The polyphase matrix associated with the filter bank in Fig. 1 is

2

E(z) = UVN −1 (z)VN −2 (z) · · · V1 (z)

(1)

where Vk (z) = IM − (1 − z −1 )vk vkT are Householder matrices; vk , 1 ≤ k < N , are unit–norm real– valued M -vectors; IM is the M × M identity matrix; and U is a M × M real–valued unitary matrix 1 2

Recall that the coding gain is the ratio of the arithmetic mean of channel variances to their geometric mean [2]. The polyphase matrix is more commonly represented as E(z) = VN −1 (z)VN −2 (z) · · · V1 (z)U [1]. While our filters

could be optimally designed based on that representation as well, the representation (1) provides a useful interpretation of U as a KLT [4], see Sec. 2.

3

[1, p. 314]. An equivalent representation of the filter bank is shown in Fig. 2 with V(z) =

Q k

Vk (z).

Let Hm (z) be the response of the mth analysis filter and H(z) be the M –vector with components {Hm (z)}. Then H(z) = E(z M )d(z) ³

where d(z) = 1, z −1 , · · · , z −M +1

´T

(2)

is a delay chain [1, p. 231].

Referring to Fig. 2, the M –vector x contains the polyphase components of the input signal x, and the M –vectors w and y are the outputs of V(z) and of the filter bank, respectively. We denote by hm the M N –vector of coefficients of Hm (z) and by um the mth column of UT .

2

The Optimization Problem

The energy compaction criterion (EC) stems from two different yet closely related formulations of the signal–adapted design problem, namely stochastic and deterministic formulations. As will be shown shortly, in each case the energy of the signal ym in the m–th channel takes the form 2 σm = hTm Rx hm ,

0≤m<M

(3)

where Rx is a certain M N × M N Toeplitz covariance matrix for the input signal x. Its first row is denoted by {rn , 0 ≤ n < M N }. Our design method would not be applicable if Rx did not have the Toeplitz property. The deterministic and stochastic formulations of the optimization problem differ 2 holds, in their definition of Rx . In each case an alternative expression for σm 2 σm = uTm Rw (0)um ,

(4)

where Rw (0) is the M × M correlation matrix for the vector w in Fig. 2.

2.1

The Stochastic Formulation

The filter bank input signal x is modeled as a wide–sense–stationary random process with Toeplitz covariance matrix 4

Rx (n, m) = E[x(n)x(m)],

0 ≤ n, m < M N.

(5)

The energy in the m–th channel is defined as the variance 4

2 σm = E|ym (n)|2

4

(6)

of the signal ym and may be written in the form (3) with Rx given by (5). Equivalently, we have 2 σm

=

Z 0.5 −0.5

|Hm (f )|2 Sx (f )df

(7)

2 is the following one. The where Sx (f ) is the spectral density of x. A third useful expression for σm

covariance of x is the M × M matrix sequence Rx (n − m) = E[x(n)xT (m)],

0 ≤ n, m < P

(8)

with components Rx,ij (n) = rM n−i+j . The covariances Rw and Ry for w and y are defined similarly. They are related to Rx by Rw (n) =

N −1 X

V(m)Rx (n + m − l)VT (l)

(9)

l,m=0

and Ry (n) = URw (n)UT

(10)

where V(m) is the (matrix–valued) impulse response of V(z). Then we obtain 2 σm = E|ym (n)|2 = Ry,mm (0) = uTm Rw (0)um ,

as announced in (4).

2.2

The Deterministic Formulation

The adaptation problem may also be formulated deterministically. Define the energy in the m–th channel as the time–averaged power of ym , 2 4 σm =

−1 1 PX |ym (p)|2 . P p=0

(11)

This definition does not require a statistical model for the signal x. We make the following two assumptions on x. (A1) The input signal x has finite length M P and is replicated by periodic extension. 0 and y (A2) The signals ym m (before and after decimation) have equal power, 1 MP

PM P −1 p=0

0 (p)|2 . |ym

5

1 P

PP −1 p=0

|ym (p)|2 =

Assumption (A2) is essentially satisfied (neglecting an O(P −1 ) stochastic term due to finite sample size) if the input signal is wide–sense stationary but is clearly much weaker. It will be approximately satisfied by any nonstationary signal with statistics varying slowly compared with the sampling interval. Under 2 still takes the form (3), with the matrix R now given by (A2), the energy σm x 4

Rx (n, m) =

P −1 1 MX x(p − n)x(p − m), M P p=0

0 ≤ n, m < M N,

namely, the empirical temporal covariance of x. By (A1), Rx is Toeplitz. Similarly, the empirical temporal covariance of the polyphase vector x is the M × M matrix sequence Rx (n − m) =

−1 1 PX x(p − n)xT (p − m), P p=0

0 ≤ n, m < P.

Again we obtain (4), with Rw (0) given by (9).

2.3

Overview of the Method

Consider the first part of (EC). The polyphase matrix E(z) is parameterized by {vm } and {um }. Referring to Fig. 2 or to (4) and (9), it is clear that the energy σ02 in the first channel depends only on {vm } and on the first row uT0 of U. A technique for maximizing σ02 is described in Sec. 2.4. This step uses (M − 1)N degrees of freedom. Once {vm } is determined, the M × M correlation matrix Rw (0) may be computed, and the only 2 , m > 0 is u . This implies that the optimal {u } are the free parameter left for maximizing σm m m

eigenvectors of the matrix Rw (0), arranged in order of decreasing eigenvalues. Thus, the second part of (EC) simply consists in performing an eigenvector decomposition of Rw (0). (Note that the first eigenvector u0 has already been identified in the first step.) As a byproduct of this operation, the optimal {um } diagonalize the output correlation matrix Ry (0) [4, Property 2]. Assume that the filter H0 (z) has maximal degree M N − 1. Then, by Theorem 4.1 in [9], the two–part technique above covers the entire set of filter banks with McMillan degree equal to M N − 1 in (2). The optimization of each Hm (f ), 1 ≤ m < M , uses M − m − 1 degrees of freedom, for a total of (M − 1)(M − 2)/2 degrees of freedom.

6

2.4

Optimization of the First Channel

The first part of (EC) is solved as follows. Define the product filter for the first channel, P0 (z) = H0 (z)H0 (z −1 ) = 1 +

MX N −1

an (z −n + z n )

(12)

n=1

where an =

X

h0 (l)h0 (l + n)

(13)

l

and aM k = 0, 1 ≤ k < N . For convenience we define the support set of the sequence (13), Λ = {n | n = kM + m, 1 ≤ m ≤ M − 1, 0 ≤ k ≤ N − 1} with cardinality |Λ| = (M − 1)N . The product filter (12) is Nyquist(M) or M th band [1], i.e., M −1 X

P0 (f + m/M ) ≡ M.

(14)

m=0

The energy in the first channel is obtained by expanding (3) and using (13) and the Toeplitz property above. The energy, normalized by 1/M for convenience 3 , is a linear function of the product filter coefficients,

"

#

X 1 2 1 E(a) = σ0 = r0 + 2 an rn . M M n∈Λ 4

(15)

It should be maximized subject to the linear constraints P0 (f ) = |H0 (f )|2 = 1 + 2

X

an cos(2πf n) ≥ 0,

0 ≤ f < 0.5.

(16)

n∈Λ

The maximization problem (15)(16) has a linear objective function with finitely many variables and infinitely many linear constraints. This is a linear SIP problem and may be thought of as an extension of the more familiar Linear Programming (LP) problems [11, 10]. Solutions to SIP present a number of desirable properties, in particular, every locally optimal solution is also globally optimal [7, 10]. The desired filter H0 (z) may be obtained from the product filter P0 (z) by spectral factorization. Each of the spectral factors of P0 (z) corresponds to a global maximum of the criterion (15). The design of the phase of H0 (z) has thus no bearing on the energy compaction properties, and we typically choose the minimum–phase solution. Also note that the matrix V(z) and the vector u0 are uniquely determined by H0 (z). The algorithm for extracting V(z) and u0 from H0 (z) follows from [9] and is summarized in Table 1. 3 1 M

PM −1 m=0

2 σm = r0 due to the paraunitary property.

7

2.5

Optimization of the Remaining Channels

Consider the second part of (EC), namely, the optimal sequential design of the filters Hm (f ), 1 ≤ m ≤ M − 1. Once the filters H0 (f ) through Hm−1 (f ) have been designed, the only free parameter left for designing Hm (f ) is um . The design of um is subject to the m + 1 orthonormality constraints uTm uk = δmk for 0 ≤ k ≤ m. There are thus only M − m − 1 degrees of freedom for designing Hm (f ). As indicated in Sec. 2.3, the optimal {um } are eigenvectors of the M ×M matrix Rw (0), arranged in order of decreasing eigenvalues. Computation of the eigenvectors is an O(M 3 ) operation. While Rw (0) may be computed using the expression (9), this operation requires N 2 M 3 multiplications, exploiting symmetries. A more efficient implementation is the following. Recursively compute the correlation matrices Rj (n) for the ouput of the Householder block Vj (z), exploiting the special structure of Vj (n), the matrix–valued impulse response of Vj (z). The recursive computation is given in Table 2 and requires only 3N (N + 1)M 2 multiplications.

3

The SIP Problem

Due to the similarity between our SIP problem (15)(16) and the SIP problem that arises in the 2–band case, we refer to [7] for more details about SIP. The problem (15)(16) is also referred to as the primal SIP, or PSIP. The set of filters a that satisfy (16) is termed the feasible set A. It inherits Properties (P1)—(P4) from [7], in particular, it is bounded, closed, convex, and symmetric around 0. We denote the maximum value of the objective function over A by V(P SIP ). The dual linear optimization problem [11] plays a central role in our work. Its variables are the Lagrange multipliers of the primal problem; once these multipliers are known, the solution may be easily determined. Additionally, the Lagrange multipliers convey useful information about the properties of the solution. The semi–infinite dual problem DSIP has as its variables generalized finite sequences λ(f ), namely the assignment of positive masses to a finite number of mass points located in the interval, 0 ≤ f ≤ 0.5. In other words, a Lagrange multiplier λ(f ) is associated with each f ; the number of nonzero multipliers is finite, but we need to consider all possible choices in the interval, 0 ≤ f ≤ 0.5. For ease of examining the duality relationships we repeat P SIP followed by its dual

8

problem. (P SIP ) maximize

E(a) =

subject to

−2

1 M

[r0 + 2

P

n∈Λ an rn ]

P

n∈Λ an cos(2πf n)

≤ 1,

0 ≤ f ≤ 0.5.

The linear semi–infinite dual problem can be stated as follows. (DSIP ) minimize

4 ¯ E(λ) =

subject to −2 and

P f

h

1 M

r0 + 2

i

P

f λ(f )

λ(f ) cos(2πnf ) = rn ,

λ(f ) ≥ 0,

n∈Λ

0 ≤ f ≤ 0.5.

We denote the minimum value of the dual problem by V(DSIP ). A basic property of DSIP is that for every a and λ satisfying the constraints of PSIP and DSIP, respectively, we have X X M ¯ [E(λ) − E(a)] = λ(f ) − an rn 2 n∈Λ f



=

X f

=

X

λ(f ) − Ã

X

X

an −2

n∈Λ

λ(f ) 1 + 2

f

X

λ(f ) cos(2πnf ) !

an cos(2πnf )

n∈Λ

f

=

 X

λ(f )P0 (f ) ≥ 0.

f

The difference V(DSIP ) − V(P SIP ), called duality gap, is thus necessarily nonnegative. The index set, 0 ≤ f ≤ 0.5, is compact, and the functions that multiply each an in (16) are continuous. This has two important implications [10]. First, approximating SIP by a finite LP through discretization of the index set [0, 0.5], will produce a solution that is arbitrarily close to an optimal SIP solution, if the discretization is made fine enough. Second, V(P SIP ) = V(DSIP ), and both PSIP and DSIP have optimal solutions that are matched through the complementary slackness condition

P f

λ(f )P0 (f ) = 0. Given an optimal solution λ(f ) to DSIP, an optimal product filter P0 (f )

may then be computed by solving a linear system of equations, as indicated in Sec. 5.

9

4

Characterization of Optimal Filters

The analysis of DSIP leads to a characterization of the optimal filters in terms of their roots on the unit circle. The fundamental result is Theorem 1 below. The proof is not given here as it parallels that of Theorem 5.1 in [7]. Condition (18) states a necessary and sufficient condition for optimality of a product filter, namely, existence of at least K ≤ (M −1)N zeroes at specified locations fk on the unit circle. The theorem does not preclude the existence of roots of P0 (z) at additional locations on the unit circle.

Theorem 1. Assume that ∃n ∈ Λ : rn 6= 0. Then there exist an optimal filter P0 (f ) with coefficients a, an integer K ∈ [1, · · · , (M − 1)N ], and certain constants ρk , fk such that ρk > 0, fk 6= fl if k 6= l, rn = −2

K X

ρk cos(2πnfk ),

n ∈ Λ,

(17)

k=1

and P0 (fk ) = 0,

1 ≤ k ≤ K.

(18)

Conditions (17) and (18) are necessary and sufficient for optimality of a ∈ A.

The K positive coefficients ρk are also the nonzero components λ(fk ) of the solution to DSIP. They can be interpreted as sensitivity parameters [7]. In particular, if the zeroes of a filter P0 (f ) are perturbed versions of the “optimal” zeroes in (18), than the objective function (15) takes the value E(a) = V(P SIP ) −

K X

λ(fk )P0 (fk ).

k=1

The solution of DSIP directly indicates which zeroes contribute significantly to coding performance. When the perturbation is small, we obtain a simple bound on the normalized error due to perturbation δfk of the optimal zeroes, V(P SIP ) − E(a) ≤ 2π 2 (M N − 1)2 max |δfk |2 . k V(P SIP ) − r0 /M

(19)

Zero Grouping. One does expect that P0 (f ) will reach its maximum possible value M for some frequencies fk∗ . Then, from the Nyquist(M) condition (14), we see that to each fk∗ must be associated M − 1 equally spaced zeroes fk∗ +

1 ∗ M , · · · , fk

+

M −1 M .

One can have up to N such frequencies fk∗ , in

which case all (M − 1)N zeroes are determined by the grouping property above. We call this property maximal zero grouping. Also interesting is the fact that the condition P0 (fk∗ ) = M implies that the 10

analysis filters H1 (f ), · · · , HM −1 (f ) necessarily have zeroes at fk∗ , owing to the Power–Complementary property

PM −1 k=0

|Hk (f )|2 = M [1, p. 290]. Our experience is that zero grouping often occurs, but we

have seen maximal zero grouping occur only when short filters are used.

Example 1. Theorem 1 has some simple but useful applications. For instance, consider the case of a sinusoidal input at frequency f ∗ , so that rn = r0 cos(2πnf ∗ ). Then any product filter satisfying the optimality conditions of Theorem 1, has M − 1 zeroes grouped at frequencies fk = f ∗ + k/M, 1 ≤ k ≤ M − 1 (hence K = M − 1), and ρk ≡ r0 /2. In this case (17) becomes the trigonometric identity M −1 X

cos(2πnm/M ) = 0,

n ∈ Λ.

k=0

These conditions are necessary and sufficient for optimality: the use of filters with additional zeroes would not improve the cost function. An equivalent condition is P0 (f ∗ ) = M . The energy in channels 1 through M is zero, due to the implied property Hm (f ∗ ) = 0 for m > 0.

Example 2. In particular, for the special case f ∗ = 0 (rn ≡ r0 , “maximally correlated” process), the necessary and sufficient condition for optimality becomes P0 (0) = M . The box filter H0 (z) =

√1 (1 M

+ z −1 + · · · + z −M +1 ) satisfies that condition.

Degeneracy. The solution to the optimization problem in each of the examples above is nonunique; the optimization problem is then said to be degenerate. Although degeneracy occurs for a very narrow class of processes, near–degeneracy frequently occurs. Examples include narrowband signals such as those considered in Sec. 7, of which Examples 1 and 2 above are limiting cases. In the 2–channel case, it has been shown that if a product filter satisfies the condition of Theorem 1 with K < (M − 1)N , then the optimization problem is degenerate, and there is a continuum of optimal filters, some of which have (M − 1)N zeroes on the unit circle [7]. We conjecture that this result holds in the M –band case.

5

SIP Algorithms

Both discretization algorithms and cutting plane algorithms [7] may be used for solving the SIP problem (15)(16). We use uniform discretization in the f –domain with P ≈ 10M N points; this is motivated by an analysis similar to that in [7, Thm 6.1] and based on the sensitivity result (19). Our 11

numerical experiments confirm that this is a reasonable design. The discretized DSIP is solved using a standard simplex algorithm. Note that Matlab’s seminf SIP optimization routine also solves linear problems using a discretization algorithm [12]. The Lagrange multipliers correspond to discrete frequencies and are typically clustered in pairs. We use a straighforward extension of Algorithm II in [7] for computing a filter P0 (z) with roots at K locations fˆk on the unit circle. This is done by assigning fˆk to the center of gravity of each pair of Lagrange multiplier indices and solving the following linear system of equations for a,   P (r) (fˆ ) = 0 k 0 (2r)  P (fˆ ) = 0 0

k

: if fˆk ∈ (0, 0.5) : else

0 ≤ r < µk , 0 ≤ k < K

where µk is the multiplicity of fˆk . The resulting product filter P0 (z) is then deflated (extracting its known K zeroes on the unit circle) so as to aid the spectral factorization routine. We used Lang and Frenzel’s spectral factorization software [13], which provides an estimate for the relative accuracy of the zeroes in the complex plane. This SIP algorithm is guaranteed to converge if the discretization is sufficiently fine.

6

Additional Design Constraints

In this section we consider modifications to our basic design so as to incorporate commonly encountered constraints. The first class of constraints deals with computational complexity. The second deals with regularity.

Computational Constraints. The implementation (1) of the filter bank is computationally attractive, since each Householder building block Vm (z) is implemented using only 2M real multiplications and one scalar delay [1, p. 315]. However, multiplication by the unitary matrix U is expensive for large M – unless U is of a special type, say a DCT matrix. The relative cost of the multiplication by U is greater for shorter filters. Because signal–adapted FIR filter banks are often capable of very high energy compaction using short filters, this raises the question, can our design be modified in a way that forces U to be a special matrix (say a DCT matrix) and optimizes energy compaction over the other degrees of freedom (Householder building blocks)? Whereas the second step of the algorithm (design of optimal columns um , 1 ≤ m ≤ M − 1) may be modified by simply using the specified um ’s, there does not seem to be a simple way to modify the first step of the algorithm which determines u0 as a 12

byproduct. We shall soon encounter a significant particular case that allows us to bypass this difficulty.

Regularity of order 1. The simplest type of regularity consists in forcing zeroes of H0 (f ) at the M − 1 frequencies f = m/M, 1 ≤ m ≤ M − 1. This implies M − 1 additional linear constraints P0 (m/M ) = 1 + 2

X

an cos(2πnm/M ) = 0,

1≤m≤M −1

(20)

n∈Λ

in the design of the product filter. These constraints are easily incorporated into the linear optimization problem (15)(16), for instance (M − 1) of the original (M − 1)N variables may be eliminated. The remaining variables are solution to a new (M − 1)(N − 1)–dimensional SIP problem. An even more convenient solution consists in retaining the original variables an and augmenting the set of linear inequalities with the M −1 equalities 2

P

n∈Λ an cos(2πnm/M )

= −1. After introducing slack variables

[11, Sec. 4.1], the dual problem becomes (DSIP 0 ) minimize

1 M

h

subject to −2 and

r0 + 2 P f

P f

λ(f ) + 2

PM −1 m=1

λ(f ) cos(2πnf ) − 2

λ(f ) ≥ 0,

i

µ(m)

PM −1 m=1

µ(m) cos(2πnm/M ) = rn ,

n∈Λ

0 ≤ f ≤ 0.5.

The additional M − 1 Lagrange multipliers µ(m) can again be interpreted as sensitivity parameters and provide information about the significance of regularity as an additional constraint. The problem DSIP is solved using a discretization algorithm similar to that in Sec. 5.

Joint Regularity and Computational Constraints. It is known that regularity of order 1 √ implies that the first row uT0 of the unitary matrix U contains constant elements 1/ M , and that the design of the remaining rows is independent of the regularity constraints [8]. But uT0 is also the first row of a Type II DCT matrix, so this property may be exploited by designing U as a Type II DCT matrix, U (k, l) = s(k) cos(πk(n + 1/2)/N ),

 p  1/M where s(k) = p  2/M

: if k = 0 : if 1 ≤ k < M,

rather than the KLT for w. In other words, if regularity and computational constraints dictate the use of a DCT matrix U, we obtain an optimal design under joint regularity and computational efficiency 13

constraints. Here the design of the first filter, H0 (f ), uses all available degrees of freedom.

Regularity of order L. This condition implies that H0 (f ) has L zeroes at frequencies f = m/M, 1 ≤ m ≤ M − 1 [8]. This implies (M − 1)L linear constraints on the product filter coefficients, µ

P

(2r)

m M

¶

X

µ

2πnm = δr0 + 2(−4π ) an n cos M n∈Λ 2 r

2r

¶

= 0,

0 ≤ r < L, 1 ≤ m < M.

(21)

Here (M − 1)L of the (M − 1)N variables a may be eliminated. The remaining variables are solution to a new (M − 1)(N − L)–dimensional SIP problem.

7

Numerical Results

The algorithm was tested on two types of input signal. In both cases, we used the stochastic formulation of the design problem in Sec. 2.1. The signals are narrowband processes, and the optimization problem may then be viewed as nearly degenerate in the sense of Sec. 4. 1. AR(1) process with correlation coefficient ρ = 0.95 (simple image model). In this case rn = ρn . 2. AR(2) process with poles at z± = ρe±iθ , with ρ = 0.975 and θ = π/3. (Models certain types of image texture.) Then rn = 2ρ cos θ rn−1 − ρ2 rn−2 , with r0 = 1 and r1 =

2ρ cos θ . 1+ρ2

The experiments reported below quantify the energy compaction performance of various filter banks. We also show frequency responses for various flavors of nonadapted and optimally adapted filters. These plots illustrate the main features that distinguish adapted and nonadapted filters, as well as the effects of regularity and/or computational constraints. In the first set of experiments, we considered the AR(1) example and compared 3–channel filter banks with length–15 filters (M = 3, N = 5), see Fig. 3 and Table 3. The nonadapted, frequency– selective filters from [1, p. 318] in Fig. 3a and the optimal adapted filters in Fig. 3b both exhibit good energy compaction performance as indicated in the corresponding columns of Table 3. The coding gain (1/M )

Q 2 1/M 2 m σm / ( m σm )

P

was 0.31 dB higher in the case of the adapted filters. Notice that

the optimal adapted filters have poor frequency selectivity, especially in the second and third channels. The optimal filters have zeroes on the unit circle, some of which are grouped as described in Sec. 4. Next, we investigated the effects of regularity constraints on the filters and compared the regular filters of Steffen et al. [8] in Fig. 3c against the optimal adapted regular filters (designed using the technique 14

in Sec. 6) in Fig. 3d. Both filter banks have good energy compaction performance. Notice from Table 3 that the regularity constraint comes at virtually no cost to energy compaction performance: the coding gain was reduced by just 0.001 dB, as the original adapted filter in Fig. 3b was already nearly regular. None of the regular filters is frequency selective. Table 4 shows the effects of filter length on energy compaction performance. The results for length– 6 (N = 2) filters are similar to those in Table 3 but somewhat worse; energy compaction properties improve with filter length. In order to assess how close the performance of the length–15 filters is to asymptotic performance, we used the expression given in [3, 4] for unconstrained–length filters m m+1 (N = ∞). These filters are brickwall filters with support set [ M , M ] for 0 ≤ m < M ; using (7) and

[14, p. 149], we obtain 2 (1/M )σm =2

Z (m+1)/M m/M

·

µ

2 1+ρ (m + 1)π mπ 1 − ρ2 df = arctan tan − tan 2 1 − 2ρ cos(2πf ) + ρ π 1−ρ 2M 2M

¶¸

These values are displayed in the last column of Table 4. The second set of experiments differs from the first in that the AR(2) example is used, see Fig. 4 and Table 5. A first, striking result is that while the energy compaction properties of the adapted filter in Fig. 4a are better than those of the nonadapted filter in Fig. 3a according to (EC), the coding gain is actually worse by 0.63 dB! Thus, the conjecture that (EC) and coding gain optimization are equivalent design criterion, is false. Also notice that while the regular nonadapted filters in Fig. 3c have poor energy compaction performance, the performance of the optimal regular adapted filters is quite reasonable: the cost of the regularity constraint is 0.18 dB in terms of coding gain. Results for length–6 filters are similar, see Fig. 5 and Table 6. In the third set of experiments, we designed 4–band filter banks adapted to the AR(2) process, see Fig. 6 and Table 7. The striking result in this case is that the frequency response of the first compaction filter H0 (f ) is not lowpass but bandpass. This is consistent with the asymptotic result (for N = ∞) which predicts that the optimal H0 (f ) is a brickwall bandpass filter with bandwidth 1/4, located around the peak of Sx (f ) at f = 1/6. The results in Table 7 show dramatic improvements over those in Table 5. In particular, the regular adapted filters outperform the regular nonadapted filters by 4.6 dB in terms of coding gain. Results for length–8 filters are shown in Table 8. In the fourth and final set of experiments, we designed 8–band filter banks adapted to the AR(2) process, see Fig. 7 and Table 9. The idea of this experiment was to investigate the performance of filters designed under computational efficiency constraints. Fig. 7a,b display the frequency response of optimal adapted and optimal adapted regular filters, respectively. Fig. 7c shows the response of 15

the filters optimized under the constraint that the matrix U in Fig. 2 is a Type II DCT matrix, using the design technique in Sec. 6. There is a 0.59 dB coding gain penalty for using regular filters, and an additional 0.39 dB penalty for using the DCT matrix instead of the optimal KLT matrix.

8

Conclusion

This study has addressed a number of basic issues in signal–adapted filter bank design using energy compaction as the adaptation criterion. Our starting point was the multistage optimization criterion (EC). The energy in the first channel depends on the coefficients of the product filter P0 (z) in (12) only, and the design of P0 (z) may be formulated as a linear SIP problem. Advantages of this formulation include: (1) every locally optimal solution is also globally optimal; (2) SIP algorithms are numerically stable; and (3) solutions to the dual problem convey useful information about statistical properties of the input signal relevant to the energy compaction problem. The particular SIP algorithm we used is based on discretization and the standard simplex method. The optimal filters H0 (z) are obtained by spectral factorization of P0 (z). While spectral factorization is a numerically ill–conditioned problem if filters are long and/or have zeroes on the unit circle, the use of deflation (removing known zeroes of P0 (z) on the unit circle) alleviates the second concern. Another advantage of signal–adapted filters is that high energy compaction performance is often obtained using short filters, so the first concern need not apply. We have been able to design length–40 filters without apparent numerical errors. In order to design the filters in the remaining channels, we used the property that the optimal unitary matrix U in the polyphase matrix representation (1) for the filter bank is a KLT matrix. The matrix V(z) in (1) is completely determined once the design of the first filter is completed. The matrix U, and hence the optimal filters in the remaining channels, are obtained using a simple matrix eigenvector decomposition. Regularity constraints on the filter H0 (z) may be formulated as additional linear constraints on the coefficients of P0 (z) and can be incorporated into the design problem. Solving this problem, one obtains globally optimal (in the energy compaction sense) filter banks with specified regularity. We have also explored the problem of optimally designing the filter bank subject to the additional constraint that the matrix U in (1) be a DCT matrix. This is simply achieved by imposing regularity of order 1 on the design of P0 (z); there are no degrees of freedom left for the design of the remaining filters. The conjecture that energy compaction and coding gain optimization are equivalent design criteria, has been shown not to hold for FIR filter banks. 16

Finally, it should be noted that the methods developed in this paper are not directly applicable to the design of signal-adapted multidimensional (M–D) FIR filter banks, due to the lack of a spectral factorization theorem for M–D polynomials. However, numerical solutions have been found for some constrained designs [15], and the solution for unconstrained–length filters is a trivial extension of known 1–D results.

17

References [1] P.P. Vaidyanathan, Multirate Systems and Filter Banks, Prentice–Hall, 1993. [2] M. Vetterli and J. Kova˘cevi´c, Wavelets and Subband Coding, Prentice-Hall, Englewood Cliffs, NJ, 1995. [3] M. K. Tsatsanis and G. B. Giannakis, “Principal Component Filter Banks for Optimal Multiresolution Analysis,” IEEE Trans. Sig. Proc., Vol. 43, No. 8, pp. 1766—1777, Aug. 1995. [4] M. Unser, “An Extension of the Karhunen-Lo`eve Transform for Wavelets and Perfect–Reconstruction Filterbanks,” SPIE Vol. 2034, pp. 45—56, 1993. [5] P. Delsarte, B. Macq and D.T.M. Slock, “Signal-Adapted Multiresolution Transform for Image Coding,” IEEE Trans. Info. Th., Vol. 38, No. 2, pp. 897–904, 1992. [6] P. Moulin, K. Ramchandran and V. Pavlovic, “Transform Image Coding based on Joint Adaptation of Filter Banks and Tree Structures,” Proc. IEEE Int. Conf. on Image Proc. (ICIP’96), Lausanne, Switzerland, Vol. 2, pp. 369—372, Sep. 1996. [7] P. Moulin, M. Anitescu, K. O. Kortanek and F. Potra, “The Role of Linear Semi–Infinite Programming in Signal–Adapted QMF Bank Design,” to appear in IEEE Trans. Sig. Proc., Sep. 1997. [8] P. Steffen, P.N. Heller, R.A. Gopinath and C.S. Burrus, “Theory of Regular M –Band Wavelet Bases,” IEEE Trans. Sig. Proc., Vol. 41, No. 12, pp. 3497—3511, 1993. [9] P.P. Vaidyanathan, T.Q. Nguyen, Z. Do˘ganata and T. Saram¨aki, “Improved Technique for Design of Perfect Reconstruction FIR QMF Banks with Lossless Polyphase Matrices,” IEEE Trans. ASSP, Vol. 37, No. 7, pp. 1042—1056, 1989. [10] R. Hettich and K.O. Kortanek, “Semi-Infinite Programming: Theory, Methods, and Applications,” SIAM Review, Vol. 35, No. 3, pp. 380–429, Sep. 1993. [11] D. G. Luenberger, Linear and Nonlinear Programming, Addison–Wesley, Reading, MA, 1984. [12] Matlab’s Optimization Toolbox, User’s Guide, The Mathworks, Inc., 1990. [13] Software available by anonymous ftp from cml.rice.edu:

c /pub/markus/software. °1992–4 LNT, Univ.

of Erlangen Nuernberg, Germany and Rice Univ., Houston, TX. [14] I.S. Gradshteyn and I.M. Ryzhik, Table of Integrals, Series, and Products, Academic Press, 1980. [15] B. Xuan and R. H. Bamberger, “Multi–Dimensional, Paraunitary Principal Component Filter Banks,” Proc. ICASSP’95, Detroit, pp. 1488–1491, 1995.

18

x(n)

- H0 (z)

y00 (n) ↓M

y0 (n) -

- H1 (z)

y10 (n) ↓M

y1 (n) -

.. .

.. .

.. .

0 -HM −1 (z) yM −1 (n)- ↓ M yM −1 (n) -

Figure 1: M –channel analysis filter bank.

x(n) ↓M

x0 (n) -

w0 (n) -

y0 (n) -

x1 (n) -

w1 (n) -

y1 (n) -

z −1

-↓ M z −1

? . ..

.. .

V(z)

.. .

U

.. .

z −1

- ↓ M xM −1 (n)-

wM −1 (n) -

yM −1 (n) -

Figure 2: Polyphase representation of filter bank in Fig. 1.

19

1

1.8

0.9

1.6

0.8

1.4

Magnitude Response

Magnitude Response

0.7 0.6 0.5 0.4

1.2 1 0.8 0.6

0.3 0.4

0.2

0.2

0.1 0 0

0.05

0.1

0.15

0.2

0.25 0.3 Frequency

0.35

0.4

0.45

0 0

0.5

0.05

0.1

0.15

0.2

1.8

1.8

1.6

1.6

1.4

1.4

1.2

1.2

1 0.8 0.6

0.2

0.15

0.2

0.25 0.3 Frequency

0.5

0.35

0.4

0.45

0.5

0.6 0.4

0.1

0.45

1

0.2

0.05

0.4

0.8

0.4

0 0

0.35

(b)

Magnitude Response

Magnitude Response

(a)

0.25 0.3 Frequency

0.35

0.4

0.45

0 0

0.5

(c)

0.05

0.1

0.15

0.2

0.25 0.3 Frequency

(d)

Figure 3: Filter response in dB for AR(1) process using 3–channel filter banks and length–15 filters (M = 3, N = 5): (a) nonadapted filters [1, p. 318], (b) adapted filters, (c) nonadapted regular filters [8], (d) adapted regular filters. Energy compaction results are shown in Table 3.

20

1.8

1.6

1.6

1.4

1.4

1.2

1.2

Magnitude Response

Magnitude Response

1.8

1 0.8 0.6

1 0.8 0.6

0.4

0.4

0.2

0.2

0 0

0.05

0.1

0.15

0.2

0.25 0.3 Frequency

0.35

0.4

0.45

0 0

0.5

0.05

0.1

0.15

(a)

0.2

0.25 0.3 Frequency

0.35

0.4

0.45

0.5

(b)

Figure 4: Filter response in dB for AR(2) process using M = 3, N = 5: (a) adapted filters, (b) adapted

1.8

1.8

1.6

1.6

1.6

1.4

1.4

1.4

1.2

1.2

1.2

1 0.8 0.6

1 0.8 0.6

0.4

0.4

0.2

0.2

0 0

0.05

0.1

0.15

0.2

0.25 0.3 Frequency

0.35

0.4

0.45

0.5

Magnitude Response

1.8

Magnitude Response

Magnitude Response

regular filters. Energy compaction results are shown in Table 5.

0 0

1 0.8 0.6 0.4 0.2

0.05

0.1

0.15

0.2

(a)

0.25 0.3 Frequency

(b)

0.35

0.4

0.45

0.5

0 0

0.05

0.1

0.15

0.2

0.25 0.3 Frequency

0.35

0.4

0.45

(c)

Figure 5: Filter response in dB for AR(2) process using M = 3, N = 2: (a) adapted filters, (b) nonadapted regular filters [8], (c) adapted regular filters. Energy compaction results are shown in Table 6.

21

0.5

2 1.8 1.6

Magnitude Response

1.4 1.2 1 0.8 0.6 0.4 0.2 0 0

0.05

0.1

0.15

0.2

0.25 0.3 Frequency

0.35

0.4

0.45

0.5

Figure 6: Filter response in dB for AR(2) process using M = 4, N = 5: adapted filters. Filter responses in successive channels are shown using solid, dotted, dashed, and dashed–dotted lines, respectively.

3

3

2.5

2.5

2.5

2

2

2

1.5

1

0.5

0 0

Magnitude Response

3

Magnitude Response

Magnitude Response

Energy compaction results are shown in Table 7.

1.5

1

0.5

0.05

0.1

0.15

0.2

0.25 0.3 Frequency

0.35

0.4

0.45

0.5

0 0

1.5

1

0.5

0.05

0.1

0.15

0.2

(a)

0.25 0.3 Frequency

(b)

0.35

0.4

0.45

0.5

0 0

0.05

0.1

0.15

0.2

0.25 0.3 Frequency

0.35

0.4

0.45

(c)

Figure 7: Filter response in dB for AR(2) process using M = 8, N = 5: (a) adapted filters, (b) adapted regular filters, (c) adapted regular/DCT filters. Energy compaction results are shown in Table 9.

22

0.5

Let CN (z) be the M –vector of polyphase components of H0 (z). Let CN (z) =

PN −1 n=0

cN (n)z −n .

for k = N − 1 downto 1 do vk =

ck+1 (k) kck+1 (k)k ,

ck (n) = ck+1 (n) − vk vkT ck+1 (n) + vk vkT ck+1 (n + 1),

0 ≤ n < k.

u0 = CN (1). Table 1: Extraction of {vk } and u0 from the first analysis filter H0 (z). Let R0 (n) = Rx (n) for 0 ≤ n ≤ N − 1. for j = 1 to N − 1 do for n = 0 to N − 1 − j do Rj (n) = Rj−1 (n) + vj vjT [Rj−1 (n + 1) − Rj−1 (n)] + [Rj−1 (n − 1) − Rj−1 (n)]vj vjT + vj vjT [2Rj−1 (n) − Rj−1 (n − 1) − Rj−1 (n + 1)]vj vjT . Rw (0) = RN −1 (0). Table 2: Recursive computation of the correlation matrix Rw (0), exploiting the special structure of the Householder blocks Vj (z). Rj (n) = RTj (−n) is the correlation matrix sequence for the output of Vj (z).

Nonadapted

Adapted

[1]

Regular

Regular

Nonadapted [8]

Adapted

σ02 /M

0.964701

0.970901

0.969977

0.970900

σ12 /M

0.025174

0.017242

0.013904

0.017213

σ22 /M

0.010125

0.011858

0.016119

0.011887

G(dB)

7.259666

7.569629

7.437999

7.568448

Table 3: Energy compaction results for AR(1) procecess using the filters shown in Fig. 3 (M = 3, N = 5).

23

N =2 Adapted

N =∞

Regular

Regular

Nonadapted [8]

Adapted

Adapted

σ02 /M

0.966843

0.966316

0.966829

0.971745

σ12 /M

0.022581

0.022296

0.022654

0.018831

σ22 /M

0.010577

0.011389

0.010517

0.009424

G(dB)

7.350620

7.262728

7.354175

7.773293

Table 4: Energy compaction results for AR(1) process using M = 3, N = 2 and M = 3, N = ∞.

Nonadapted

Adapted

[1]

Regular

Regular

Nonadapted [8]

Adapted

σ02 /M

0.510620

0.511788

0.487307

0.511730

σ12 /M

0.481925

0.475783

0.225747

0.474162

σ22 /M

0.007454

0.012428

0.286946

0.014108

G(dB)

4.350480

3.625736

0.231372

3.447308

Table 5: Energy compaction results for AR(2) process using the filters shown in Fig. 4 (M = 3, N = 5).

Adapted

Regular

Regular

Nonadapted [8]

Adapted

σ02 /M

0.507783

0.470706

0.507074

σ12 /M

0.487144

0.469587

0.487425

σ22 /M

0.005073

0.059707

0.005501

G(dB)

4.900063

1.493825

4.784172

Table 6: Energy compaction results for AR(2) process using the filters shown in Fig. 5 (M = 3, N = 2).

24

Adapted

Regular

Regular

Nonadapted [8]

Adapted

σ02 /M

0.946137

0.080306

0.834606

σ12 /M

0.034523

0.071440

0.150786

σ22 /M

0.014504

0.781062

0.011744

σ32 /M

0.004836

0.067192

0.002864

G(dB)

8.079426

2.782665

7.412881

Table 7: Energy compaction results for AR(2) process using M = 4, N = 5. The adapted filters are shown in Fig. 6.

Adapted

Regular

Regular

Nonadapted [8]

Adapted

σ02 /M

0.918218

0.151022

0.809100

σ12 /M

0.065887

0.639997

0.133339

σ22 /M

0.012225

0.187483

0.053898

σ32 /M

0.003670

0.021498

0.003664

G(dB)

7.895271

2.502952

5.658318

Table 8: Energy compaction results for AR(2) process using M = 4, N = 2. The adapted filters are shown in Fig. 6.

25

Adapted

Regular

Regular&DCT

Adapted

Adapted

σ02 /M

0.898201

0.880607

0.880493

σ12 /M

0.048480

0.046782

0.024187

σ22 /M

0.026532

0.033070

0.029521

σ32 /M

0.011648

0.021391

0.039093

σ42 /M

0.008946

0.011939

0.019361

σ52 /M

0.003268

0.003227

0.004017

σ62 /M

0.001742

0.001799

0.002108

σ72 /M

0.001183

0.001186

0.001220

G(dB)

9.832958

9.245101

8.854587

Table 9: Energy compaction results for AR(2) process using M = 8, N = 5.

26

List of Figures 1

M –channel analysis filter bank. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

2

Polyphase representation of filter bank in Fig. ??. . . . . . . . . . . . . . . . . . . . . .

19

3

Filter response in dB for AR(1) process using 3–channel filter banks and length–15 filters (M = 3, N = 5): (a) nonadapted filters [1, p. 318], (b) adapted filters, (c) nonadapted regular filters [8], (d) adapted regular filters. Energy compaction results are shown in Table ??. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

Filter response in dB for AR(2) process using M = 3, N = 5: (a) adapted filters, (b) adapted regular filters. Energy compaction results are shown in Table ??. . . . . . . .

5

20

21

Filter response in dB for AR(2) process using M = 3, N = 2: (a) adapted filters, (b) nonadapted regular filters [8], (c) adapted regular filters. Energy compaction results are shown in Table ??. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

21

Filter response in dB for AR(2) process using M = 4, N = 5: adapted filters. Filter responses in successive channels are shown using solid, dotted, dashed, and dashed– dotted lines, respectively. Energy compaction results are shown in Table ??. . . . . . .

7

22

Filter response in dB for AR(2) process using M = 8, N = 5: (a) adapted filters, (b) adapted regular filters, (c) adapted regular/DCT filters. Energy compaction results are shown in Table ??. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

22

List of Tables 1

Extraction of {vk } and u0 from the first analysis filter H0 (z). . . . . . . . . . . . . . .

2

Recursive computation of the correlation matrix Rw (0), exploiting the special structure

23

of the Householder blocks Vj (z). Rj (n) = RTj (−n) is the correlation matrix sequence for the output of Vj (z). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

23

Energy compaction results for AR(1) procecess using the filters shown in Fig. ?? (M = 3, N = 5). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

4

Energy compaction results for AR(1) process using M = 3, N = 2 and M = 3, N = ∞.

24

5

Energy compaction results for AR(2) process using the filters shown in Fig. ?? (M = 3, N = 5). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

Energy compaction results for AR(2) process using the filters shown in Fig. ?? (M = 3, N = 2). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

9

24

Energy compaction results for AR(2) process using M = 4, N = 5. The adapted filters are shown in Fig. ??. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

24

25

Energy compaction results for AR(2) process using M = 4, N = 2. The adapted filters are shown in Fig. ??. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25

Energy compaction results for AR(2) process using M = 8, N = 5. . . . . . . . . . . .

26

28