Nonuniform Multirate Filter Banks: Analysis and Design ... - CiteSeerX
Nonuniform Multirate Filter Banks: Analysis and Design with an H1 Performance Measure Tongwen Chen Dept. of Electrical and Computer Engineering University of Calgary Calgary, Alberta, Canada T2N 1N4 Email: [email protected] July 31, 1995 Abstract
A general multirate lter-bank system with multiple channels and nonuniform bands is studied with the reconstruction performance measured by the worst-case energy gain (H1 norm) of the error system between the multirate system and a pure time-delay system. Using blocking and polyphase decomposition, we associate with the multirate system an equivalent linear time-invariant (LTI) system whose transfer matrix can be computed by a simple procedure. Based on this LTI system, the optimal design problem can be reduced to one of H1 optimization, which can be solved by ready-made software in many cases. For illustration, optimal synthesis lters are designed for a two-channel example and their properties are studied in detail.
EDICS Number: SP 2.4.1. Keywords: digital signal processing, multirate ltering, nonuniform lter banks, fractional decimation, H1 optimization.
1 Introduction This paper treats a general multirate system depicted in Figure 1. Shown are m channels numbered 0; 1; ; m ? 1, each channel consisting of an analysis subsystem and a synthesis subsystem. The i-th analysis subsystem has a fractional decimation ratio of pi =qi , where pi and qi are coprime integers, and is realized by the the pi -fold expander " pi (upsampling by pi ), a linear-time invariant (LTI) band-pass lter Hi , and the qi -fold decimator # qi (downsampling by qi ). In a typical application, e.g., subband coding, the incoming discrete-time signal x is split into m di erent frequency bands by the analysis subsystems and then the subband signals
This research was supported by the Natural Sciences and Engineering Research Council of Canada.
1
x-
0
- "p ? r 1
rr
0
- #q
0
v0-
"q
0
-
F0
- #p
H1
- #q
1
v1-
"q
1
-
F1
- # p - ?j ?r
"p - H
- " pm? - Hm? - # qm? 1
1
1
vm?-1
0
1
" qm? - Fm? - # pm? 1
1
1
rr - ?jy-
Figure 1: The general multirate lter system.
vi are coded and transmitted; ignoring the coding and transmission errors, we arrive at the synthesis end, whose purpose is to reconstruct x as close as possible from the subband signals. The i-th synthesis subsystem converts the decimated signal back to the normal frequency and has a fractional decimation ratio of qi =pi ; this is the combination of " qi , Fi , and # pi , where Fi is an LTI lter, suitably designed. The outputs of the synthesis subsystems are summed up to form the reconstructed signal y . For ease of reference, the integers qi and pi are called the decimation and expansion inte-
gers, respectively, with reference to the analysis end. Throughout the paper, we assume the analysis bank fHig and the synthesis bank fFi g consist of LTI systems with real-coecient transfer functions. The setup in Figure 1 allows arbitrary (fractional) nonuniform band splitting to be accomplished; for example, by suitable choice of the analysis lters, the frequency band [0; ) can be split into the following m subbands:
"
!
mX ?2 p mX ?1 p i; i : 0; pq 0 ; pq 0 ; pq 0 + pq 1 ; ; q q 0 0 0 1 i=0 i i=0 i
(1)
For the bands not to overlap, we assume mX ?1 p i i=0 qi
= 1:
(2)
This also preserves the sampling density. Thus, the i-th channel processes the frequency components of x in the i-th band of (1). By proper choice of integers pi and qi , one can practically match any requirement for band splitting. 2
Nonuniform lter banks are especially useful in digital audio [5, 27]. Perfect-reconstruction nonuniform lter banks were studied before: Hoang and Vaidyanathan [9] treated a nonuniform system and gave a necessary condition for cases in which alias cancellation is impossible to achieve; Nayebi et. al. [20] developed a time-domain design approach for FIR lters; Kovacevic and Vetterli [14] accomplished design of nonuniform lter banks by transforming them into uniform systems; nally, Liu and Bruton [16] treated perfect reconstruction in a somewhat di erent setup involving certain modulators. Note that the m-channel uniform lter-bank system is a special case of the setup in Figure 1 by setting all pi to 1 and all qi to m; in this case the literature is quite extensive, see the book by Vaidyanathan [27] and the references therein. In this paper we take a di erent approach to nonuniform lter bank design. Instead of designing for perfect reconstruction, we design for close-to-perfect reconstruction. This idea was explored in the uniform lter banks [13, 15, 21, 4]. From the study in [4], relaxing the condition for perfect reconstruction allows freedom to use any reasonable analysis lters to satisfy the band-splitting and coding requirements; then the synthesis lters are designed to get as close to perfect reconstruction as possible. Under some mild condition, one can always get arbitrarily close to perfect reconstruction by trading o lter complexity and time delay in reconstruction. To quantify the degree of closeness to perfect reconstruction, we use the worst-case energy gain as in [4] for the uniform lter banks. For this we need the Hilbert space `2 of squaresummable discrete-time signals, real and perhaps vector-valued, de ned over the time set of nonnegative integers. The inner product and norm on `2 are
hu; vi =
X
k
u(k)0v(k);
kuk = 2
X
k
ku(k)k
2
!1=2
;
where 0 denotes transpose and the norm on u(k) is the Euclidean one. For perfect reconstruction to occur in Figure 1, we must have y (k) = x(k ? d) for some delay integer d 0. Thus for close-to-perfect reconstruction systems, we should compare y (k) with x(k ? d) and de ne the reconstruction error signal
e(k) := x(k ? d) ? y(k): This error of course depends on the input x in Figure 1. To eliminate this dependence, de ne the reconstruction performance measure in Figure 1 to be
J := sup kkxekk2 ; 2
3
where the supremum is taken over all nonzero x in `2 . Note that J = 0 implies perfect reconstruction and that if J is small, the reconstruction error, measured by `2 norm, is uniformly small over all possible inputs x with unit norm. J is an H1 performance index for reasons to be seen in Section 4.2. Let the system x 7! y in Figure 1 be T . For xed delay integer d, de ne the ideal delay system to be Tid with transfer function z ?d . Then the quantity J is exactly the `2 -induced norm of the operator Tid ? T , J = kTid ? T k: (3) J is shown to be a good indicator for the reconstruction performance of lter banks [4, 18]. The focus of this paper is to answer the following analysis and design problems related to the nonuniform lter system in Figure 1:
Given the system in Figure 1 and the associated delay integer d, how to compute J in general? (analysis)
Given the decimation and expansion integers pi and qi, the analysis lters Hi, and the delay integer d, how to design stable, causal synthesis lters Fi to minimize J ? (design) These two problems were dealt with for uniform lter banks (with two channels) in [4] by the frequency-domain H1 model-matching theory [7, 8, 3]. To extend the results to nonuniform banks, the major diculty lies in that there is no frequency-domain model available for the general system in Figure 1. Developing such a model is then the rst task of this paper. Note that the design approach takes the viewpoint that the analysis lters are xed and the synthesis lters are to be designed to minimize J . This viewpoint are useful in the following scenarios:
Analysis lters have previously been designed by coding engineers to satisfy their coding and band splitting requirements.
Analysis lters are designed based on other methods and then use H1 optimization to redesign the synthesis lters to improve reconstruction performance (J ). Fixing analysis lters, design synthesis lters to optimize J ; then x the synthesis lters and redesign the analysis lters to optimize J ; and repeat several iterations. The advantage is that every optimization is guaranteed solvable and the global optimality attainable.
The paper is organized as follows. Section 2 presents the blocking technique used to transform linear periodic systems into equivalent LTI ones as well as its frequency-domain 4
interpretation. Section 3 derives via blocking the equivalent LTI system for Figure 1; the transfer matrix of this LTI system can be constructed simply using polyphase decompositions of the analysis and synthesis lters. Based on the framework established in Section 3, Section 4 presents some analysis results for the multirate system in Figure 1 and discusses a solution to the analysis problem stated earlier. In Section 5, the design problem is reduced to an H1 -optimal model matching problem, which can be solved directly in many cases by known techniques in control; moreover, design examples are also discussed in detail for illustration. Finally, Section 6 points out a design constraint called structural dependency and o ers some concluding remarks. (A preliminary conference version of this paper appeared in [2].) We conclude this section by introducing some notation. Both time-varying and timeinvariant systems occur in this paper; they are represented as time-domain operators, denoted by capital letters, e.g., G and T . If a system G is LTI, its transfer function (matrix) is written G^(z). For ease of use in equations, we also write the p-fold decimator (# p) as Dp and the q-fold expander (" q) as Eq. Thus the system x 7! y in Figure 1 can be simply written
T=
mX ?1 i=0
(Dp Fi Eq ) (Dq Hi Ep ) : i
i
i
i
This is the sum of m subsystems, each channel representing a subsystem. The i-th subsystem is in fact periodic with period qi , or simply, qi -periodic. This implies that the overall system is n-periodic, where n is the least common multiple of all qi , i.e.,
n = l.c.m.fq0; q1; ; qm?1 g:
(4)
This fact suggests the technique of blocking, which was used for lter banks with uniform bands [28] and with a generalized block sampling scheme [11].
2 Blocking Signals and Systems The standard technique for treating periodic/multirate systems is called blocking in signal processing [17, 19, 28] and lifting in control [12]. This section collects some facts about blocking and relates it to polyphase decomposition in the frequency domain.
2.1 Blocking Signals Let ` be the space of discrete-time signals, real and perhaps vector valued, de ned on the time set f0; 1; 2; g. A signal v in ` can be written
v = fv(0); v(1); v(2); g: 5
For an integer n > 0, de ne the n-fold blocking operator, Ln , via v = Ln v (underlining denotes blocking): 9 82 3 3 2 > > v ( n ) v (0) > > > > 7 6 > > 7 = 6 .. 7 ; 6 (5) .. 7 ; > : > > . . 5 5 4 4 > > > > ; : v (n ? 1) v(2n ? 1) Ln maps ` to `n , the external direct sum of n copies of `. If the underlying period for v is h, the underlying period for the blocked signal v is nh. Note that the blocking operation results in no loss of information. The inverse L?n 1 , mapping `n to `, amounts to reversing the operation in (5). In the frequency domain, bring in the n-fold polyphase decomposition for the z -transform of v : nX ?1 v^(z) = z?j v^j (z n): j =0
This representation is unique for a given integer n. Now de ne the n-fold polyphase vector of v^ as 2 3 v^0(z) 6 v^1(z) 777 6 6 .. 7: 6 . 5 4 v^n?1 (z) Summarized below are some properties of blocking, which are straightforward to show.
Lemma 1 Let v = Ln v. (a) The z -transform of v , v^(z ), equals the n-fold polyphase vector of v^(z ). (b) Ln is norm-preserving on `2 , i.e., kvk2 = kv k2.
Lemma 1 (a) characterizes the blocking operation in the frequency domain.
2.2 Blocking Systems For a linear discrete-time system G mapping ` to `, blocking of its input and output signals induces another linear system, denoted G, as is depicted in the commutative diagram in Figure 2. The system G is the blocked system and we get from the diagram G = Ln GL?n 1 . If G maps v to y and is single-input and single-output (SISO), G maps v to y and has n inputs and n outputs. It is clear that there exists a one-to-one correspondence between G and G.
Lemma 2 Assume G is a linear SISO system with input v and output y. Let G = LnGL?n . 1
6
G
`n Ln
6 G
`
-
`n
-
`
6Ln
Figure 2: The commutative diagram of blocking. (a) G is time-invariant i G is n-periodic. (b) If G is n-periodic, the transfer matrix G^ of G is n n and relates to the input and output polyphase vectors v^ and y^ by y^(z ) = G^ (z )^v(z ). (c) Blocking preserves the `2 -induced norm, i.e., kGk = kGk.
The results in this lemma are standard [17, 28, 12]. Part (a) is the basis for using blocking in periodic systems; part (b) suggests a simple way to compute frequency responses of periodic systems; part (c) will be useful for handling our performance index J , an `2-induced norm. In Section 3.4, we shall see that a periodic system can be implemented by its blocked, LTI system using decimators and expanders. If G is already LTI, of course G is LTI. How do we write down the transfer matrix of G based on the transfer function of G?
Lemma 3 The following two statements are equivalent: (a) G is LTI with the transfer function
G^ (z) =
nX ?1 z?j G^ j (z n): j =0
(6)
(b) G is LTI with the transfer matrix 2 6 6 6 ^ G(z) = 666 6 4
G^ 0(z) z ?1 G^ n?1 (z) G^ 1(z) G^ 0(z) G^ 2(z) G^ 1 (z) .. .
G^n?1 (z)
.. .
G^ n?2 (z) 7
z ?1 G^ n?2 (z) z ?1 G^ n?1 (z) G^ 0(z) .. .
G^ n?3 (z)
3
z ?1 G^ 1(z) 7 z ?1 G^ 2(z) 77 z ?1 G^ 3(z) 77 : 7 .. 7 . 5 G^ 0(z)
(7)
One use of this lemma is to test if a periodic system reduces to an LTI system, and hence has the alias-free property, by looking at the blocked system. The direction (a) ) (b) is well-known and can be found in, e.g., [12]; the converse can be proven easily. For an LTI system G, the representation in (6) is the n-fold polyphase decomposition of G; the n functions G^j are the n-fold polyphase components of G. The matrix in (7) is referred to as the n-fold polyphase matrix associated with G. Note that this n n matrix is Toeplitz, namely, constant along all diagonals, and is determined by only n functions. Another version of the result in Lemma 3 is characterized by a property called pseudocirculant in [28]. A pseudocirculant matrix does not have exactly the same form in (7); this is because the blocking operation in [28] is slightly di erent from the one de ned in (5): The blocked vectors contain the same set of input values but have di erent permutations.
3 Blocking Fractionally Decimated Filter Banks The overall system T (x 7! y ) in Figure 1 is the sum of m analysis-synthesis subsystems,
T=
mX ?1 i=0
Si Ai ;
(8)
where Ai and Si are the i-th analysis and synthesis subsystems:
Ai = Dq Hi Ep ; Si = Dp Fi Eq : i
i
i
i
(9)
As we commented before, the i-th subsystem Si Ai is qi -periodic and hence T is n-periodic, where n is given in (4). In this section we wish to block T to get a frequency-domain model; this requires blocking all subsystems Si Ai by the n-fold blocking operator.
3.1 Blocking Analysis Subsystems First, let us consider how to block appropriately the analysis subsystems Ai . Since all Ai have a similar structure, in this subsection we drop the subscripts and de ne A = Dq HEp, as depicted in Figure 3. The results in this subsection are of independent interest too, because
x-
Ep
-
H
-
Dq
v-
Figure 3: An analysis subsystem. 8
the setup in Figure 3 represents a general system which changes sampling rates by a factor of p=q . The system A : x 7! v in Figure 3 inherits some properties from the original setup in Figure 1: The system H is LTI; and we wish to do n-fold blocking to the input x, where n is a multiple of q. How to block the output v? It is not hard to see that in general A is not n-periodic due to the fact that A changes sampling periods from, say, h, to qh=p. For the blocked system to have the same rate at the input and output, we should block v by the (np=q )-fold blocking operator. (Note that np=q is an integer.) This suggests blocking A this way, A := Lnp=q ALn?1 : Indeed, this A is LTI. To see this, let U be the unit time-delay system on ` and U the unit time-advance system (the transfer functions for U and U are z ?1 and z , respectively). To show that A is LTI, it suces to show U AU = A, i.e., A is shift-invariant. By the obvious identities, U Dq = Dq (U )q and Ep U = U p Ep , and time invariance of H , we get (U )np=q AU n = = = =
(U )np=q Dq HEpU n Dq (U )np HU npEp
Dq HEp A:
This fact together with the identities U Lnp=q = Lnp=q (U )np=q and L?n 1 U = U n L?n 1 further implies
U AU = U Lnp=q AL?n 1 U = Lnp=q (U )np=q AU n L?n 1 = Lnp=q AL?n 1 = A; and hence A is LTI and has a transfer matrix. Now de ne two linear transformations Jp and Kq as block diagonal matrices: 82 3 2 3 9 > > 1 1 > > > > > 7 6 7 > = 66 .. 77 ; ; 66 .. 77 > ; > 4 . 5 > 4 . 5 > > > > ; : 0 0 p1 p1
Kq = block diag
h
i
h
i
1 0 0 1q ; ; 1 0 0 1q : 9
The blocks in Jp are repeated n times and the blocks in Kq (np=q ) times. Thus Jp is (np) n and Kq is (np=q ) (np). The following lemma is useful in obtaining the transfer matrix of A.
Lemma 4 The following identities hold: Lnp=q Dq = Kq Lnp; EpL?n 1 = L?np1 Jp : This means that when combined with blocking operations properly, decimation and expansion amount to multiplication by certain matrices. The lemma can be veri ed easily from the de nitions of the quantities involved. Based on Lemma 4, we have
A = Lnp=q Dq HEpL?n 1 = Kq (Lnp HL?np1)Jp : The system Lnp HL?np1 is the (np)-fold blocking of H and so the transfer matrix for A is ^ p; A^ = Kq HJ where H^ , with size (np) (np), is the (np)-fold polyphase matrix associated with H and is given in Lemma 3. Note that pre- and post-multiplying H^ by Kq and Jp , respectively, amount to extracting an (np=q ) n matrix from H^ by taking rows numbered 0; q; 2q; ; np ? q and columns numbered 0; p; 2p; ; np ? p. The result is summarized below.
Theorem 1 Let H^ be the (np)-fold polyphase matrix of H : 3 2 ^ H (z) z ? H^ np? (z) z ? H^ (z) 6 H^ (z) H^ (z) z? H^ (z) 777 6 ^ 6 0
H (z) = 6 4
1
.. .
H^ np?1 (z)
1
1
1
1
0
.. .
H^ np?2 (z)
1
.. 7: 5 . ^ H0(z) 2
Then the (np=q ) n transfer matrix for A is obtained by retaining every q -th row and every p-th column of H^ starting from the upper-left corner: 3 2 ^ H0(z) z ?1 H^ np?p (z) z ?1 H^ p(z) 6 ^ (z ) H^ q?p (z) z ?1 H^ q+p (z) 777 6 H A^(z) = 66 q.. .. .. 7 5 4 . . . H^ np?q (z) H^ np?q?p (z) H^ np?q+p (z) (assuming q > p).
10
Note that the matrix A^ is no longer Toeplitz. The frequency-domain interpretation of A^ can be obtained by Lemma 1 as follows: v^(z ) = A^(z )^x(z ), where x^ and v^ are the n-fold and (np=q )-fold polyphase vectors of x^(z ) and v^(z ) in Figure 3, respectively. In general, n is a multiple of q . If n = q , to get the matrix A^(z ), we need to perform the (pq )-fold polyphase decomposition of H ; such a fact was also observed by Hsiao [10] from a somewhat di erent point of view. Can we recover H^ (z ) from A^(z )? Since blocking is one-to-one, the equivalent question is: Can we recover H from Dq HEp?
Lemma 5 Assume H is LTI. Then the mapping H 7! Dq HEp is one-to-one if p and q are coprime.
The converse of the lemma is also true, i.e., if p and q are not coprime, then H 7! Dq HEp is not one-to-one.
Proof of Lemma 5 It suces to show that Dq HEp = 0 implies H = 0. So assume
Dq HEp = 0. Let h be the impulse response of H , i.e., h = H, where is the unit impulse signal. It follows that for any integer i, Dq HEpU i = 0, where U is the unit time-delay system. This implies, by the de nition of Ep , Dq HU ip = 0, and, by time invariance of H and the de nition of Dq , h(ip + jq) = 0; 8i; j: (10) Since p and q are coprime, there exist integers p and q such that qp + pq = 1. Thus, for any k, take i = kq and j = kp in (10) to get h(k) = 0; and hence H = 0. 2
Based on Lemma 5, we conclude that the matrix A^(z ) in Theorem 1 contains all (np)-fold polyphase components of H if p and q are coprime.
3.2 Blocking Synthesis Subsystems This subsection relates to blocking the synthesis subsystems Si . Let S be the dual of the system A studied in Section 3.1, namely, the system v 7! y shown in Figure 4, with F LTI and the same assumptions about integers p; q; n: n q > p and n is a multiple of q . We wish to block y by Ln . Since S changes the sampling rates by a factor q=p, we de ne the blocked S to be 1 S = LnSL?np=q : Similarly as in Section 3.1, we can show that S is LTI and the transfer matrix is given by ^ q; S^ = KpFJ 11
v-
-
Eq
-
F
Dp
y-
Figure 4: A synthesis subsystem. where F^ is the (np)-fold polyphase matrix of F and the matrices Jq and Kp are de ned by 82 > > > > 66 > 4 > > : h
2
3
3
9
> 1 1 > > 7 6 > = 0 0 777 7 6 7 6 ; Jq = block .. 7 ; ; 6 .. 7 > 4 . 5 > .5 > > 0 q1 ; 0 q1 h i i Kp = block diag 1 0 0 1p ; ; 1 0 0 1p : Jq has (np=q) diagonal blocks and the size is (np) (np=q); Kp has n diagonal blocks and the size is n (np). In this way we arrive at a dual result.
Theorem 2 Let F^ be the (np)-fold polyphase matrix of F : 2 ^ 3 F (z) z ? F^np? (z) z ? F^ (z) 6 F^ (z) F^ (z) z? F^ (z) 777 6 ^ 6 0
F (z) = 6 4
1
.. .
F^np?1 (z)
1
1
1
1
0
.. .
F^np?2 (z)
1
.. 7: 5 . ^ F0 (z) 2
Then the n (np=q ) transfer matrix for S is obtained by retaining every p-th row and every q-th column of F^ starting from the upper-left corner: 2 ^ 3 F0(z) z ?1 F^np?q (z) z ?1 F^q (z) 6 ^ 7 ?1 ^ ?1 ^ ^S (z ) = 666 Fp.(z ) z Fnp.?q+p (z ) z Fp.+q (z ) 777 .. .. .. 4 5 ^ ^ ^ Fnp?p (z) Fnp?p?q (z) Fq?p (z) (assuming q > p).
Again, in general S^ is no longer Toeplitz but it contains all (np)-fold polyphase components of F if p and q are coprime (Lemma 5). Similar frequency-domain interpretation of S^ can be made by Lemma 1. 12
3.3 Blocking the Overall System Now we shall put together the results in the preceding two subsections to get the n-fold blocking of the overall system T . First, block the analysis and synthesis subsystems Ai and Si in (9) and de ne
Ai = Lnp =q Ai L?n 1 ; Si = LnSi L?np1 =q ; i = 0; 1; ; m ? 1: i
i
i
i
By Theorems 1 and 2, Ai and Si are LTI and their transfer matrices can be computed, respectively, from the (npi )-fold polyphase decompositions of Hi and Fi : A^i is (npi =qi ) n and S^i is n (npi=qi ). Assuming the two matrices A^i and S^i are computed, we proceed to compute T = Ln TL?n 1 . By (8),
T = = =
mX ?1 i=0 mX ?1 i=0 mX ?1 i=0
Ln SiAi L?n 1 (Ln SiL?np1 =q )(Lnp =q Ai L?n 1 ) i
i
i
i
Si Ai:
Hence the transfer matrix for T is
T^(z) =
mX ?1 i=0
S^i (z)A^i (z):
This corresponds to the blocked system in Figure 5, where all systems involved are LTI and signals are blocked, e.g., vi is the (npi =qi )-fold blocking of vi in Figure 1. De ne the analysis matrix and the synthesis matrix to be, respectively, 2 ^ 3 A0(z) 6 ^ 7 h i A^(z) = 664 A1(z) 775 ; S^(z) = S^0(z) S^1 (z) S^m?1 (z) : A^m?1 (z) Note that the analysis (synthesis) matrix consists of polyphase components of the analysis (synthesis) bank only. Thus we conclude
T^(z) = S^(z)A^(z): 13
x
? r
rr
-
A0
-
A1
v0
-
S0
v1
-
S1
vm?1
- Am?
1
- ?j ?r
rr - ?j y -
- Sm?
1
Figure 5: The blocked LTI system. Both the analysis and synthesis matrices, A^ and S^, are n n; for example, since A^i is (npi =qi ) n, the number of columns of A^ is n and the number of rows is
n
mX ?1 p i i=0 qi
= n;
by the assumption in (2). Therefore, the blocked T has a transfer matrix which can be expressed as the product of the synthesis and analysis matrices. In the frequency domain, this means 2 6 4
y^0(z) .. .
y^n?1 (z)
3 2 7 = S^(z )A^(z )6 5 4
x^0 (z) .. .
x^n?1 (z)
3 7; 5
where x^j and y^j are the n-fold polyphase components of the input x and output y in Figure 1 by Lemma 1. Now let us summarize the procedure to compute the n n matrix T^ in Figure 5:
Step 1 For i = 0; 1; ; m ? 1, compute the (npi)-fold polyphase decompositions of Hi and Fi :
H^ i (z) =
npX i ?1 z?j H^ ij (z npi ); j =0
F^i (z) =
npX i ?1 z?j F^ij (znpi ): j =0
Step 2 For i = 0; 1; ; m ? 1, form the (npi)-fold polyphase matrices, H^ i and F^i, associated with Hi and Fi , respectively (Lemma 3); then de ne A^i to be the (npi =qi ) n 14
matrix obtained by retaining every qi -th row and every pi-th column of H^ i , starting from the upper-left corner (Theorem 1), 2 3 H^ i0(z) z ?1 H^ i;np ?p (z) z ?1 H^ ip (z) 6 H^ iq (z) H^ i;q ?p (z) z ?1 H^ i;q +p (z) 777 6 ^ 6 Ai (z) = 6 (11) .. .. .. 7; 4 5 . . . H^ i;np ?q (z) H^ i;np ?q ?p (z) H^ i;np ?q +p (z) i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
and similarly, de ne S^i to be the n (npi =qi ) matrix obtained by retaining every pi -th row and every qi -th column of F^i , starting from the upper-left corner (Theorem 2), 2 6 ^Si (z ) = 666 4
F^i0 (z) F^ip (z)
z ?1 F^i;np ?q (z) z ?1 F^iq (z) ? z 1 F^i;np ?q +p (z) z ?1 F^i;p +q (z) i
.. .
i
F^i;np ?p (z) i
.. .
i
i
i
i
i
F^i;np ?p ?q (z) i
i
.. .
i
i
i
i
F^i;q ?p (z) i
3 7 7 7 7: 5
(12)
i
Step 3 Form the two n n analysis and synthesis matrices A^(z) = S^(z) = and then
h h
i0 A^00(z) A^1 0 (z) A^m?1 0(z) ; i S^0(z) S^1 (z) S^m?1 (z)
T^(z) = S^(z)A^(z):
The transfer matrix T^(z ) serves as a frequency-domain model for the general multirate system in Figure 1. The analysis and synthesis matrices A^(z ) and S^(z ) can be regarded as generalization of the polyphase representation matrices for uniform lter banks [29, 26]. Since we assumed that pi and qi are coprime, there exist one-to-one correspondences between A and fHig and between S and fFi g (Lemma 5). In [14] nonuniform lter banks are studied by converting them into uniform lter banks via block diagram manipulations. The frequency-domain model derived here amounts to accomplish the same but uses only matrix manipulations. One can perhaps argue that the method in this paper is more systematical and concise. However, the time-domain approach in [20] is quite di erent. Finally, we remark that a di erent approach for frequency-domain modeling of multirate systems is based on alias-component matrices; this approach was developed for certain multirate systems involving fractional sampling factors in [25, 24]; see also [23]. 15
3.4 Implementation via Blocking We shall conclude Section 3 by looking at an alternative structure, based on the n n analysis and synthesis matrices, to implement the multirate system T in Figure 1. Recall that the cascaded LTI system S A equals the blocked T , namely, Ln TL?n 1 ; hence T = L?n 1 S ALn . It turns out that Ln and L?n 1 can be realized by the n-fold decimator and expander together with the time delay U and advance U .
Lemma 6 The following identities hold: 2 6 6 Ln = 66 4
Dn Dn U .. .
Dn (U )n?1
3 7 7 7 7; 5
i
h
L?n 1 = En UEn U n?1 En :
(13)
The proof is straightforward from the de nitions of the quantities involved. Based on this lemma, we have h
T = En UEn
2 6 i 6 n ? 1 U En S A66 4
Dn DnU .. .
Dn(U )n?1
3 7 7 7 7: 5
This corresponds to the system shown in Figure 6. Note that in this structure, the z transforms of the signals x0 ; x1; ; xn?1 are the n-fold polyphase components of the input x in Figure 1 and similarly for y0 ; y1; ; yn?1 and the output y; the LTI systems A and S operate at a rate which is slower than the normal rate by a factor of n. Such a structure can be used to implement any periodic system based on its blocked, LTI system.
4 Analysis In this section we deal with some analysis issues associated with the multirate system in Figure 1; in particular, we give a solution to the analysis problem stated in Section 1.
4.1 Conditions for Alias Free and Perfect Reconstruction The overall system T in Figure 1 in general su ers aliasing, magnitude distortion, and phase distortion. The system is alias-free if T is LTI and has perfect reconstruction if, in addition, T is free from magnitude and phase distortions; in this case, T equals some pure time-delay 16
x
?
- #n
x0-
y0-
-
"n
z
? r
- #n
x1-
rr
rr r
A
?
y1-
-
z?1
"n
xn?-1
- 6j 6r
rr
S
6
z
- #n
-j y6
yn?1
-
z?1
"n
6
Figure 6: The implementation structure via blocking. system. The conditions for alias free and perfect reconstruction can be stated in terms of the analysis and synthesis matrices, A and S , derived in the preceding section.
Theorem 3 Let T be the nonuniform lter-bank system in Figure 1. (a) T is alias-free i there exists some functions G^ j (z ), j = 0; 1; ; n ? 1, such that 2 ^ 3 G z? G^n? z ? G^ 6 ^ G^ z? G^ 777 ^ ^ 66 G 0
S A = 6 .. 4 . G^ n?1
1
1
1
1
.. .
1
0
G^ n?2
1
.. 7 : . 5 ^ G0 2
(14)
Moreover, if this condition is satis ed, the transfer function for T is
T^(z) =
nX ?1 z?j G^ j (z n ): j =0
(b) T achieves perfect reconstruction i there exists some integers k and r, k 0, 0 r n ? 1, such that " # ?1 Ir 0 z ? k S^ A^ = z ; (15) I 0 n?r
where Ir and In?r are the r r and (n ? r) (n ? r) identity matrices, respectively. Moreover, if this condition is satis ed, T^(z ) = z ?(kn+r) .
17
Proof Part (a) follows directly from Lemma 3. For part (b), recall that T has perfect
reconstruction i T = G, where G is some pure time-delay system: G^ (z ) = z ?d for some d 0. This condition is equivalent to Ln TL?n 1 = LnGL?n 1 , or in the frequency domain, S^(z)A^(z) = G^ (z). To get G^(z), the n-fold polyphase matrix for G, write d = kn + r for some unique integers k and r with k 0 and 0 r n ? 1; then the polyphase components of G are given by ( ?k j = r; ^ Gj (z) = 0z j 6= r: Hence G^ (z ) equals to the right-hand side of (15) by Lemma 3 and the proof is complete. 2 The conditions given in this theorem have been obtained in other forms in [20, 14].
4.2
H1
Analysis
Now we answer the analysis question raised in Section 1 regarding computation of the performance measure J . This is done by evaluating the H1 norm of some transfer matrix. To review, the Hardy space H1 is the Banach space of complex-valued matrix functions of z that are analytic and bounded outside the unit disk, that is, for jz j > 1. Thus H1 is the space of transfer matrices of stable, causal, LTI discrete-time systems. The norm of a matrix G^(z) in H1 is its supremal maximum singular value on the unit circle: h i kG^k1 = sup max G^ ej! : !
If G^ (z ) is the transfer matrix of a stable, causal, LTI discrete-time system G, then the `2 induced norm of G equals the H1 norm of G^ (z ), that is, kGk = kG^ k1 . Back to our problem of computing J . Given the nonuniform system T in Figure 1 and the ideal delay system Tid with T^id (z ) = z ?d for some d 0, the problem is to compute the `2-induced norm in (3), i.e, J = kTid ? T k. This cannot be evaluated directly since T is not LTI. However, blocking Tid ? T preserves the induced norm [Lemma 2 (c)]. So
J = kTid ? T k:
(16)
Now the systems Tid and T are both LTI and so the induced norm in (16) equals the H1 norm of the transfer matrix, J = kT^id ? S^ A^k1; where T^id is derived in Theorem 3 (b): "
0
T^id(z) = z ?k I n?r
#
z ?1 Ir ; d = kn + r; 0 r n ? 1: 0 18
(17)
In summary, given a nonuniform system, computing the reconstruction performance index J is reduced to evaluating the H1 norm of some matrix function involving the analysis and synthesis matrices. The latter problem is standard [7, 8, 4] and MATLAB software exists [1].
5 Design This section studies the optimal design problem stated in Section 1.
5.1
H1
-Optimal Design
For the optimal design problem associated with Figure 1 we assume the analysis lters Hi have been designed based on band-splitting and coding requirements and now we are to design the synthesis lters Fi to minimize J in (3), where T^id (z ) = z ?d and the delay integer d is assumed to be given. Since blocking preserves the induced norm (as in Section 4.2), design of Fi relates to design of the synthesis matrix S ; the equivalent design problem based on the frequency-domain model is now: Given T^id and A, design S^ to minimize J = kT^id ? S^ A^k1 , the H1 norm of the n n transfer matrix T^id ? S^ A^. The matrix T^id depends on the delay integer d and is given in (17); the analysis matrix A^ depends on Hi and is constructed by the procedure in Section 3.3. If S^ can be designed, the synthesis lters Fi are recovered according to the same procedure. The equivalent problem is an H1 optimization problem:
Jopt = inf^ kT^id ? S^ A^k1 : S
(18)
For xed analysis systems, the optimal performance Jopt in general is a function of the delay integer d; so we may write Jopt (d). Computing this function allows us to study the tradeo between reconstruction performance and reconstruction delay. It is a fact proven in [4] for uniform lter banks that lim J (d) = 0: d!1 opt In systems with the same decimation integer, namely, with all equal qi in Figure 1, there is no structural dependency in S^, see Section 6 for a discussion. In this case, we optimize over the class of S^ in (18) which are just stable, causal, and IIR; this is a standard H1 -optimal model-matching problem, which can be solved by standard techniques in control [7, 8, 3] and there exists MATLAB software [1]. For illustration, consider the example below. 19
5.2 Design Example Shown in Figure 7 is a two-channel nonuniform system, in which q0 = q1 = 3, p0 = 2, and
x-
"2 - H
0
- #3
v0-
"3 - F
-
H1
- #3
v1-
"3 - F
- #2
0
? y-j
1
Figure 7: The two-channel example.
p1 = 1. The frequency range [0; ) is split into the low-frequency band [0; 2=3) for channel 0 and the high-frequency band [2=3; ) for channel 1. To get the frequency-domain model, we note that n = 3 and then follow the procedure in Section 3.3. For the analysis matrix, bring in the 6-fold and 3-fold polyphase decompositions of H0 and H1 , respectively,
H^ 0(z) = to get
5 X z?j H^ 0j (z6); j =0
2 A^(z) = 64
H^ 1(z) =
2 X z?j H^ 1j (z3); j =0
3
H^ 00(z) H^ 03(z) H^ 10(z)
z ?1 H^ 04(z) z ?1 H^ 02(z) H^ 01(z) z ?1 H^ 05(z) 75 : z ?1 H^ 12(z) z ?1 H^ 11(z) Similarly, the polyphase decompositions of F0 and F1 , X X F^0(z) = z ?j F^0j (z6); F^1 (z) = z ?j F^0j (z3);
yield the synthesis matrix
5
2
j =0
j =0
2 S^(z) = 64
F^00(z) F^02(z) F^04(z)
z ?1 F^03(z) z ?1 F^05(z) F^01(z)
(19)
3
F^10(z) F^11(z) 75 : F^12(z)
(20)
For an ideal delay system T^id (z ) = z ?d , represent the blocked system Tid in the form of (17) with n = 3. Then the associated H1 optimization problem is Jopt(d) = inf kT^id ? S^ A^k1 : S^
20
If d, H0, and H1 are given, so are T^id and A^; and then the optimization problem can be solved by using the MATLAB function hinfsyn in the -Analysis and Synthesis Toolbox [1]. In this case, the discrete-time problem needs to be converted into an equivalent continuous-time problem via bilinear transformation, see, e.g., [3], and all computations involved are based on state-space data. In this way, the optimal S computed and therefore the optimal F0 and F1 obtained from (20) and (19) are in general IIR. For illustration, we compute the function Jopt (d) for two sets of FIR linear-phase analysis lters, H0 lowpass with cuto frequency =3 (this agrees with the fact that channel 0 processes frequencies up to 2=3 because the 2-fold expander is prior to H0 ) and H1 highpass with cuto frequency 2=3, all obtained by the MATLAB function r1 which uses the Hamming window: Case 1: order of H0 = 19 and order of H1 = 20; Case 2: order of H0 = 23 and order of H1 = 24. The curves for Jopt (d) are shown in Figure 8. The tradeo between reconstruction performance 0
10
−1
10
−2
10
−3
10
−4
10
0
5
10
15
20
25
30
Figure 8: Jopt versus d for Case 1 (solid) and Case 2 (dash). and delay is obvious: For small Jopt , one needs to tolerate relatively large d. It can be proven as in [4] that lim d!1 Jopt (d) = 0 under a mild regularity condition on the analysis matrix, which means that for reasonable analysis lters, one can get arbitrarily close to perfect reconstruction (Jopt = 0) by tolerating suciently large time delay. 21
For a closer look at the optimal synthesis lters computed, we consider the analysis lters in Case 1. The magnitude Bode plots for H0 and H1 are given in Figure 9. Select d = 20. 20
0
−20
−40
−60
−80
−100
−120 0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
Figure 9: Case 1: jH^ 0j (solid) and jH^ 1j (dash) in dB versus !=2 . From the solid curve in Figure 8 we nd Jopt = 0:8066%. The corresponding optimal synthesis matrix S^(z ) is IIR with order 40, from which we compute from (20) and (19) the IIR synthesis lters with their impulse responses shown in Figure 10. The IIR impulse responses, despite their di erent e ective lengths, look very symmetric, indicating that the designed F0 and F1 yield almost linear phase. This suggests to approximate them simply by truncation. Now we truncate the impulse responses, f0 (k) and f1 (k), of F0 and F1 to get their FIR approximations: f0FIR (k) = f0(k); 5 k 56; f1FIR (k) = f1 (k); 0 k 23: With the approximated lters F0FIR and F1FIR , there is a price to be paid in reconstruction performance. Using the analysis result in Section 4.2, the performance achieved by F0FIR and F1FIR is computed to be JFIR = 1:132% (compare with Jopt = 0:8066%). This means that for all input x, the energy gain from x to the reconstruction error e, kek22 =kxk22, is always less than or equal to (1:132%)2 = 0:000128. Figure 11 gives the Bode plots for F0FIR (lowpass) and F1FIR (highpass). With the approximated FIR synthesis lters, the reconstruction errors, e(k) := x(k ? 20) ? y (k), are simulated for x(k) = 1(k), the unit step signal, and x(k) = cos(0:5k) + cos(2:5k), the sum of a low-frequency sinusoid (! = 0:5) and a high-frequency sinusoid (! = 2:5), see Figure 12 for the error signals. 22
3 2 1 0 −1 0
10
20
30
40
50
60
10
20
30
40
50
60
1.5 1 0.5 0 −0.5 −1 0
Figure 10: Impulse responses f0 (k) and f1 (k) versus k for the optimal F0 (upper) and F1 (lower).
20 10 0 −10 −20 −30 −40 −50 −60 −70 −80 0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
Figure 11: jF^0FIRj (solid) and jF^1FIR j (dash) in dB versus !=2 . 23
0.02
0.015
0.01
0.005
0
−0.005
−0.01
−0.015
−0.02 0
10
20
30
40
50
60
70
80
90
100
Figure 12: e(k) versus k for the unit step input (solid) and the sinusoid input (dot). Finally, we remark that better reconstruction performance is possible by using larger time delays, longer synthesis lters, or even other analysis lters.
6 Conclusion The method used in the preceding section for IIR H1 optimal lter design extends directly to the class of nonuniform lter systems with equal decimation integers, namely, the class of systems in Figure 1 with all qi equal. However, if qi are not all the same, there is a diculty. Consider the three-channel case with p0 ; p1 ; p2 = 1 ; 1 ; 1 : q0 q1 q2 2 4 4 In this case n = 4. To write down the synthesis matrix, we introduce the 4-fold polyphase decompositions for Fi , 3 X F^i (z) = z?j F^ij (z 4 ); i = 0; 1; 2; j =0
then the procedure in Section 3.3 yields 2 ^ F00 z?1 F^02 6 ?1 ^ ^ S^ = 664 FF^01 z F^F03 02 00 ^ ^ F03 F01 24
F^10 F^11 F^12 F^13
F^20 3 F^21 777 : F^22 5 F^23
Note that the rst two columns in the synthesis matrix S^ are related to each other: Each is obtained from the four polyphase components of F0 . Such a phenomenon is called structural dependency, which poses a constraint in design. The design approach in Section 5 cannot handle this constraint easily because there the H1 -optimal synthesis matrices are computed from explicit formulas based on solving two Riccati equations [6]. It can be shown that structural dependency occurs whenever qi are not all equal, under the assumption in (2). We remark that this phenomenon was also observed in [14]. One way to handle structural dependency in H1 design is to use numerical optimization. In [22], convex optimization was applied to design H1 -optimal FIR lters in QMF banks. This approach is promising for treating structural dependency in nonuniform lter bank design.
Acknowledgement:
This paper was revised and expanded based on an earlier manuscript submitted to IEEE Trans. on Signal Processing in 1994. The author would like to thank the Associate Editor, Professor Roberto Bamberger, and an anonymous reviewer for encouraging this revision, and to the reviewers for helpful comments.
References [1] G. J. Balas, J. C. Doyle, K. Glover, A. Packard, and R. Smith, -Analysis and Synthesis Toolbox User's Guide, MuSyn Inc. and the MathWorks, Inc., 1991. [2] T. Chen, \A lifting approach to analysis and design of nonuniform multirate lter banks," Proc. Canadian Conf. on Electrical and Computer Engineering, Montreal, Sept. 1995, pp. 874-877. [3] T. Chen and B. A. Francis, Optimal Sampled-Data Control Systems, Springer, London, 1995. [4] T. Chen and B. A. Francis, \Design of multirate lter banks by H1 optimization," IEEE Trans. Signal Processing, vol. 43, pp. 2822-2830, 1995. [5] R. E. Crochiere and L. R. Rabiner, Multirate Digital Signal Processing, Prentice-Hall, Englewood Cli s, NJ, 1983. [6] J. C. Doyle, K. Glover, P. P. Khargonekar, and B. A. Francis, \State-space solutions to standard H2 and H1 problems," IEEE Trans. Automatic Control, vol. 34, no. 8, pp. 831-847, 1989. [7] B. A. Francis, A Course in H1 Control Theory, Springer-Verlag, Lecture Notes in Control and Information Sciences, vol. 88, Berlin, 1987. [8] M. Green and D. J. N. Limebeer, Linear Robust Control, Prentice-Hall, Englewood Cli s, NJ, 1995. 25
[9] P.-Q. Hoang and P. P. Vaidyanathan, \Non-uniform multirate lter banks: theory and design," Proc. IEEE Int. Symp. Circuits and Syst., Portland, OR, 1989, pp. 371-374. [10] C.-C. Hsiao, \Polyphase lter matrix for rational sampling rate conversions," Proc. IEEE Int. Conf. Accoust., Speech, Signal Processing, Dallas, April 1987, pp. 2173-2176. [11] M. R. K. Khansari and A. Leon-Garcia, \Subband decomposition of signals with generalized sampling," IEEE Trans. Signal Processing, vol. 41, no. 12, pp. 3365-3376, 1993. [12] P. P. Khargonekar, K. Poolla, and A. Tannenbaum, \Robust control of linear timeinvariant plants using periodic compensation," IEEE Trans. Automatic Control, vol. 30, pp. 1088-1096, 1985. [13] R. D. Koilpillai and P. P. Vaidyanathan, \A spectral factorization approach to pseudoQMF design," IEEE Trans. Signal Processing, vol. 41, no. 1, pp. 82-92, 1993. [14] J. Kovacevic and M. Vetterli, \Perfect reconstruction lter banks with rational sampling factors," IEEE Trans. Signal Processing, vol. 41, no. 6, pp. 2047-2064, 1993. [15] Y. C. Lim, R. H. Yang, and S.-N. Koh, \The design of weighted minimax quadrature mirror lters," IEEE Trans. Signal Processing, vol. 41, no. 5, pp. 1780-1789, 1993. [16] B. Liu and L. T. Bruton, \The design of nonuniform-band maximally decimated lter banks," Proc. IEEE Int. Symp. Circuits and Syst., Chicago, IL, 1993, pp. 375-378. [17] R. A. Meyer and C. S. Burrus, \A uni ed analysis of multirate and periodically timevarying digital lters," IEEE Trans. Circuits and Systems, vol. 22, pp. 162-168, 1975. [18] S. Mirabbasi, B. A. Francis, and T. Chen, \Controlling distortions in maximally decimated lter banks," IEEE Trans. on Circuits and Systems II: Analog and Digital Signal Processing, to appear. [19] T. Miyawaki and C. W. Barnes, \Multirate recursive digital lters { a general approach and block structures," IEEE Trans. Acoustics, Speech, Signal Processing, vol. 31, pp. 1148-1154, 1983. [20] K. Nayebi, T. P. Barnwell III, and M. J. T. Smith, \Nonuniform lter banks: a reconstruction and design theory," IEEE Trans. Signal Processing, vol. 41, no. 3, pp. 1114-1127, 1993. [21] T. Q. Nguyen, \Near-perfect-reconstruction pseudo-QMF banks," IEEE Trans. Signal Processing, vol. 42, no. 1, pp. 65-76, 1994. [22] S. Ratzla , An Optimal Model Matching Approach to Multirate FIR Filter Bank Design, M.Sc. Thesis, Dept. of Electrical and Computer Engineering, University of Calgary, August 1995. [23] R. L. Reng, \Polyphase and modulation descriptions of multirate systems { a systematic approach," Proc. Int. Conf. DSP, 1995, pp. 212-217. 26
[24] R. G. Shenoy, \Analysis of multirate components and application to multirate lter design," Proc. IEEE Int. Conf. Accoust., Speech, Signal Processing, 1994, vol. 3, pp. 121-124. [25] R. G. Shenoy, D. Burnside, and T. W. Parks, \Linear periodic systems and multirate lter design," IEEE Trans. Signal Processing, vol. 42, pp. 2242-2256, 1994. [26] P. P. Vaidyanathan, \Theory and design of M -channel maximally decimated Quadrature mirror lters with arbitrary M , having the perfect-reconstruction property," IEEE Trans. Acoustics, Speech, Signal Processing, vol. 35, no. 4, pp. 476-492, 1987. [27] P. P. Vaidyanathan, Multirate Systems and Filter Banks, Prentice-Hall, Englewood Cli s, NJ, 1993. [28] P. P. Vaidyanathan and S. K. Mitra, \Polyphase networks, block digital ltering, LPTV systems, and alias-free QMF banks: a uni ed approach based on pseudocirculants," IEEE Trans. Acoustics, Speech, Signal Processing, vol. 36, no. 3, pp. 381-391, 1988. [29] M. Vetterli, \A theory of multirate lter banks," IEEE Trans. Acoustics, Speech, Signal Processing, vol. 35, no. 3, pp. 356-372, 1987.
27
A general multirate lter-bank system with multiple channels and nonuniform bands is studied with the reconstruction performance measured by the worst-case energy gain (H1 norm) of the error system between the multirate system and a pure time-delay system. Using blocking and polyphase decomposition, we associate with the multirate system an equivalent linear time-invariant (LTI) system whose transfer matrix can be computed by a simple procedure. Based on this LTI system, the optimal design problem can be reduced to one of H1 optimization, which can be solved by ready-made software in many cases. For illustration, optimal synthesis lters are designed for a two-channel example and their properties are studied in detail.
EDICS Number: SP 2.4.1. Keywords: digital signal processing, multirate ltering, nonuniform lter banks, fractional decimation, H1 optimization.
1 Introduction This paper treats a general multirate system depicted in Figure 1. Shown are m channels numbered 0; 1; ; m ? 1, each channel consisting of an analysis subsystem and a synthesis subsystem. The i-th analysis subsystem has a fractional decimation ratio of pi =qi , where pi and qi are coprime integers, and is realized by the the pi -fold expander " pi (upsampling by pi ), a linear-time invariant (LTI) band-pass lter Hi , and the qi -fold decimator # qi (downsampling by qi ). In a typical application, e.g., subband coding, the incoming discrete-time signal x is split into m di erent frequency bands by the analysis subsystems and then the subband signals
This research was supported by the Natural Sciences and Engineering Research Council of Canada.
1
x-
0
- "p ? r 1
rr
0
- #q
0
v0-
"q
0
-
F0
- #p
H1
- #q
1
v1-
"q
1
-
F1
- # p - ?j ?r
"p - H
- " pm? - Hm? - # qm? 1
1
1
vm?-1
0
1
" qm? - Fm? - # pm? 1
1
1
rr - ?jy-
Figure 1: The general multirate lter system.
vi are coded and transmitted; ignoring the coding and transmission errors, we arrive at the synthesis end, whose purpose is to reconstruct x as close as possible from the subband signals. The i-th synthesis subsystem converts the decimated signal back to the normal frequency and has a fractional decimation ratio of qi =pi ; this is the combination of " qi , Fi , and # pi , where Fi is an LTI lter, suitably designed. The outputs of the synthesis subsystems are summed up to form the reconstructed signal y . For ease of reference, the integers qi and pi are called the decimation and expansion inte-
gers, respectively, with reference to the analysis end. Throughout the paper, we assume the analysis bank fHig and the synthesis bank fFi g consist of LTI systems with real-coecient transfer functions. The setup in Figure 1 allows arbitrary (fractional) nonuniform band splitting to be accomplished; for example, by suitable choice of the analysis lters, the frequency band [0; ) can be split into the following m subbands:
"
!
mX ?2 p mX ?1 p i; i : 0; pq 0 ; pq 0 ; pq 0 + pq 1 ; ; q q 0 0 0 1 i=0 i i=0 i
(1)
For the bands not to overlap, we assume mX ?1 p i i=0 qi
= 1:
(2)
This also preserves the sampling density. Thus, the i-th channel processes the frequency components of x in the i-th band of (1). By proper choice of integers pi and qi , one can practically match any requirement for band splitting. 2
Nonuniform lter banks are especially useful in digital audio [5, 27]. Perfect-reconstruction nonuniform lter banks were studied before: Hoang and Vaidyanathan [9] treated a nonuniform system and gave a necessary condition for cases in which alias cancellation is impossible to achieve; Nayebi et. al. [20] developed a time-domain design approach for FIR lters; Kovacevic and Vetterli [14] accomplished design of nonuniform lter banks by transforming them into uniform systems; nally, Liu and Bruton [16] treated perfect reconstruction in a somewhat di erent setup involving certain modulators. Note that the m-channel uniform lter-bank system is a special case of the setup in Figure 1 by setting all pi to 1 and all qi to m; in this case the literature is quite extensive, see the book by Vaidyanathan [27] and the references therein. In this paper we take a di erent approach to nonuniform lter bank design. Instead of designing for perfect reconstruction, we design for close-to-perfect reconstruction. This idea was explored in the uniform lter banks [13, 15, 21, 4]. From the study in [4], relaxing the condition for perfect reconstruction allows freedom to use any reasonable analysis lters to satisfy the band-splitting and coding requirements; then the synthesis lters are designed to get as close to perfect reconstruction as possible. Under some mild condition, one can always get arbitrarily close to perfect reconstruction by trading o lter complexity and time delay in reconstruction. To quantify the degree of closeness to perfect reconstruction, we use the worst-case energy gain as in [4] for the uniform lter banks. For this we need the Hilbert space `2 of squaresummable discrete-time signals, real and perhaps vector-valued, de ned over the time set of nonnegative integers. The inner product and norm on `2 are
hu; vi =
X
k
u(k)0v(k);
kuk = 2
X
k
ku(k)k
2
!1=2
;
where 0 denotes transpose and the norm on u(k) is the Euclidean one. For perfect reconstruction to occur in Figure 1, we must have y (k) = x(k ? d) for some delay integer d 0. Thus for close-to-perfect reconstruction systems, we should compare y (k) with x(k ? d) and de ne the reconstruction error signal
e(k) := x(k ? d) ? y(k): This error of course depends on the input x in Figure 1. To eliminate this dependence, de ne the reconstruction performance measure in Figure 1 to be
J := sup kkxekk2 ; 2
3
where the supremum is taken over all nonzero x in `2 . Note that J = 0 implies perfect reconstruction and that if J is small, the reconstruction error, measured by `2 norm, is uniformly small over all possible inputs x with unit norm. J is an H1 performance index for reasons to be seen in Section 4.2. Let the system x 7! y in Figure 1 be T . For xed delay integer d, de ne the ideal delay system to be Tid with transfer function z ?d . Then the quantity J is exactly the `2 -induced norm of the operator Tid ? T , J = kTid ? T k: (3) J is shown to be a good indicator for the reconstruction performance of lter banks [4, 18]. The focus of this paper is to answer the following analysis and design problems related to the nonuniform lter system in Figure 1:
Given the system in Figure 1 and the associated delay integer d, how to compute J in general? (analysis)
Given the decimation and expansion integers pi and qi, the analysis lters Hi, and the delay integer d, how to design stable, causal synthesis lters Fi to minimize J ? (design) These two problems were dealt with for uniform lter banks (with two channels) in [4] by the frequency-domain H1 model-matching theory [7, 8, 3]. To extend the results to nonuniform banks, the major diculty lies in that there is no frequency-domain model available for the general system in Figure 1. Developing such a model is then the rst task of this paper. Note that the design approach takes the viewpoint that the analysis lters are xed and the synthesis lters are to be designed to minimize J . This viewpoint are useful in the following scenarios:
Analysis lters have previously been designed by coding engineers to satisfy their coding and band splitting requirements.
Analysis lters are designed based on other methods and then use H1 optimization to redesign the synthesis lters to improve reconstruction performance (J ). Fixing analysis lters, design synthesis lters to optimize J ; then x the synthesis lters and redesign the analysis lters to optimize J ; and repeat several iterations. The advantage is that every optimization is guaranteed solvable and the global optimality attainable.
The paper is organized as follows. Section 2 presents the blocking technique used to transform linear periodic systems into equivalent LTI ones as well as its frequency-domain 4
interpretation. Section 3 derives via blocking the equivalent LTI system for Figure 1; the transfer matrix of this LTI system can be constructed simply using polyphase decompositions of the analysis and synthesis lters. Based on the framework established in Section 3, Section 4 presents some analysis results for the multirate system in Figure 1 and discusses a solution to the analysis problem stated earlier. In Section 5, the design problem is reduced to an H1 -optimal model matching problem, which can be solved directly in many cases by known techniques in control; moreover, design examples are also discussed in detail for illustration. Finally, Section 6 points out a design constraint called structural dependency and o ers some concluding remarks. (A preliminary conference version of this paper appeared in [2].) We conclude this section by introducing some notation. Both time-varying and timeinvariant systems occur in this paper; they are represented as time-domain operators, denoted by capital letters, e.g., G and T . If a system G is LTI, its transfer function (matrix) is written G^(z). For ease of use in equations, we also write the p-fold decimator (# p) as Dp and the q-fold expander (" q) as Eq. Thus the system x 7! y in Figure 1 can be simply written
T=
mX ?1 i=0
(Dp Fi Eq ) (Dq Hi Ep ) : i
i
i
i
This is the sum of m subsystems, each channel representing a subsystem. The i-th subsystem is in fact periodic with period qi , or simply, qi -periodic. This implies that the overall system is n-periodic, where n is the least common multiple of all qi , i.e.,
n = l.c.m.fq0; q1; ; qm?1 g:
(4)
This fact suggests the technique of blocking, which was used for lter banks with uniform bands [28] and with a generalized block sampling scheme [11].
2 Blocking Signals and Systems The standard technique for treating periodic/multirate systems is called blocking in signal processing [17, 19, 28] and lifting in control [12]. This section collects some facts about blocking and relates it to polyphase decomposition in the frequency domain.
2.1 Blocking Signals Let ` be the space of discrete-time signals, real and perhaps vector valued, de ned on the time set f0; 1; 2; g. A signal v in ` can be written
v = fv(0); v(1); v(2); g: 5
For an integer n > 0, de ne the n-fold blocking operator, Ln , via v = Ln v (underlining denotes blocking): 9 82 3 3 2 > > v ( n ) v (0) > > > > 7 6 > > 7 = 6 .. 7 ; 6 (5) .. 7 ; > : > > . . 5 5 4 4 > > > > ; : v (n ? 1) v(2n ? 1) Ln maps ` to `n , the external direct sum of n copies of `. If the underlying period for v is h, the underlying period for the blocked signal v is nh. Note that the blocking operation results in no loss of information. The inverse L?n 1 , mapping `n to `, amounts to reversing the operation in (5). In the frequency domain, bring in the n-fold polyphase decomposition for the z -transform of v : nX ?1 v^(z) = z?j v^j (z n): j =0
This representation is unique for a given integer n. Now de ne the n-fold polyphase vector of v^ as 2 3 v^0(z) 6 v^1(z) 777 6 6 .. 7: 6 . 5 4 v^n?1 (z) Summarized below are some properties of blocking, which are straightforward to show.
Lemma 1 Let v = Ln v. (a) The z -transform of v , v^(z ), equals the n-fold polyphase vector of v^(z ). (b) Ln is norm-preserving on `2 , i.e., kvk2 = kv k2.
Lemma 1 (a) characterizes the blocking operation in the frequency domain.
2.2 Blocking Systems For a linear discrete-time system G mapping ` to `, blocking of its input and output signals induces another linear system, denoted G, as is depicted in the commutative diagram in Figure 2. The system G is the blocked system and we get from the diagram G = Ln GL?n 1 . If G maps v to y and is single-input and single-output (SISO), G maps v to y and has n inputs and n outputs. It is clear that there exists a one-to-one correspondence between G and G.
Lemma 2 Assume G is a linear SISO system with input v and output y. Let G = LnGL?n . 1
6
G
`n Ln
6 G
`
-
`n
-
`
6Ln
Figure 2: The commutative diagram of blocking. (a) G is time-invariant i G is n-periodic. (b) If G is n-periodic, the transfer matrix G^ of G is n n and relates to the input and output polyphase vectors v^ and y^ by y^(z ) = G^ (z )^v(z ). (c) Blocking preserves the `2 -induced norm, i.e., kGk = kGk.
The results in this lemma are standard [17, 28, 12]. Part (a) is the basis for using blocking in periodic systems; part (b) suggests a simple way to compute frequency responses of periodic systems; part (c) will be useful for handling our performance index J , an `2-induced norm. In Section 3.4, we shall see that a periodic system can be implemented by its blocked, LTI system using decimators and expanders. If G is already LTI, of course G is LTI. How do we write down the transfer matrix of G based on the transfer function of G?
Lemma 3 The following two statements are equivalent: (a) G is LTI with the transfer function
G^ (z) =
nX ?1 z?j G^ j (z n): j =0
(6)
(b) G is LTI with the transfer matrix 2 6 6 6 ^ G(z) = 666 6 4
G^ 0(z) z ?1 G^ n?1 (z) G^ 1(z) G^ 0(z) G^ 2(z) G^ 1 (z) .. .
G^n?1 (z)
.. .
G^ n?2 (z) 7
z ?1 G^ n?2 (z) z ?1 G^ n?1 (z) G^ 0(z) .. .
G^ n?3 (z)
3
z ?1 G^ 1(z) 7 z ?1 G^ 2(z) 77 z ?1 G^ 3(z) 77 : 7 .. 7 . 5 G^ 0(z)
(7)
One use of this lemma is to test if a periodic system reduces to an LTI system, and hence has the alias-free property, by looking at the blocked system. The direction (a) ) (b) is well-known and can be found in, e.g., [12]; the converse can be proven easily. For an LTI system G, the representation in (6) is the n-fold polyphase decomposition of G; the n functions G^j are the n-fold polyphase components of G. The matrix in (7) is referred to as the n-fold polyphase matrix associated with G. Note that this n n matrix is Toeplitz, namely, constant along all diagonals, and is determined by only n functions. Another version of the result in Lemma 3 is characterized by a property called pseudocirculant in [28]. A pseudocirculant matrix does not have exactly the same form in (7); this is because the blocking operation in [28] is slightly di erent from the one de ned in (5): The blocked vectors contain the same set of input values but have di erent permutations.
3 Blocking Fractionally Decimated Filter Banks The overall system T (x 7! y ) in Figure 1 is the sum of m analysis-synthesis subsystems,
T=
mX ?1 i=0
Si Ai ;
(8)
where Ai and Si are the i-th analysis and synthesis subsystems:
Ai = Dq Hi Ep ; Si = Dp Fi Eq : i
i
i
i
(9)
As we commented before, the i-th subsystem Si Ai is qi -periodic and hence T is n-periodic, where n is given in (4). In this section we wish to block T to get a frequency-domain model; this requires blocking all subsystems Si Ai by the n-fold blocking operator.
3.1 Blocking Analysis Subsystems First, let us consider how to block appropriately the analysis subsystems Ai . Since all Ai have a similar structure, in this subsection we drop the subscripts and de ne A = Dq HEp, as depicted in Figure 3. The results in this subsection are of independent interest too, because
x-
Ep
-
H
-
Dq
v-
Figure 3: An analysis subsystem. 8
the setup in Figure 3 represents a general system which changes sampling rates by a factor of p=q . The system A : x 7! v in Figure 3 inherits some properties from the original setup in Figure 1: The system H is LTI; and we wish to do n-fold blocking to the input x, where n is a multiple of q. How to block the output v? It is not hard to see that in general A is not n-periodic due to the fact that A changes sampling periods from, say, h, to qh=p. For the blocked system to have the same rate at the input and output, we should block v by the (np=q )-fold blocking operator. (Note that np=q is an integer.) This suggests blocking A this way, A := Lnp=q ALn?1 : Indeed, this A is LTI. To see this, let U be the unit time-delay system on ` and U the unit time-advance system (the transfer functions for U and U are z ?1 and z , respectively). To show that A is LTI, it suces to show U AU = A, i.e., A is shift-invariant. By the obvious identities, U Dq = Dq (U )q and Ep U = U p Ep , and time invariance of H , we get (U )np=q AU n = = = =
(U )np=q Dq HEpU n Dq (U )np HU npEp
Dq HEp A:
This fact together with the identities U Lnp=q = Lnp=q (U )np=q and L?n 1 U = U n L?n 1 further implies
U AU = U Lnp=q AL?n 1 U = Lnp=q (U )np=q AU n L?n 1 = Lnp=q AL?n 1 = A; and hence A is LTI and has a transfer matrix. Now de ne two linear transformations Jp and Kq as block diagonal matrices: 82 3 2 3 9 > > 1 1 > > > > > 7 6 7 > = 66 .. 77 ; ; 66 .. 77 > ; > 4 . 5 > 4 . 5 > > > > ; : 0 0 p1 p1
Kq = block diag
h
i
h
i
1 0 0 1q ; ; 1 0 0 1q : 9
The blocks in Jp are repeated n times and the blocks in Kq (np=q ) times. Thus Jp is (np) n and Kq is (np=q ) (np). The following lemma is useful in obtaining the transfer matrix of A.
Lemma 4 The following identities hold: Lnp=q Dq = Kq Lnp; EpL?n 1 = L?np1 Jp : This means that when combined with blocking operations properly, decimation and expansion amount to multiplication by certain matrices. The lemma can be veri ed easily from the de nitions of the quantities involved. Based on Lemma 4, we have
A = Lnp=q Dq HEpL?n 1 = Kq (Lnp HL?np1)Jp : The system Lnp HL?np1 is the (np)-fold blocking of H and so the transfer matrix for A is ^ p; A^ = Kq HJ where H^ , with size (np) (np), is the (np)-fold polyphase matrix associated with H and is given in Lemma 3. Note that pre- and post-multiplying H^ by Kq and Jp , respectively, amount to extracting an (np=q ) n matrix from H^ by taking rows numbered 0; q; 2q; ; np ? q and columns numbered 0; p; 2p; ; np ? p. The result is summarized below.
Theorem 1 Let H^ be the (np)-fold polyphase matrix of H : 3 2 ^ H (z) z ? H^ np? (z) z ? H^ (z) 6 H^ (z) H^ (z) z? H^ (z) 777 6 ^ 6 0
H (z) = 6 4
1
.. .
H^ np?1 (z)
1
1
1
1
0
.. .
H^ np?2 (z)
1
.. 7: 5 . ^ H0(z) 2
Then the (np=q ) n transfer matrix for A is obtained by retaining every q -th row and every p-th column of H^ starting from the upper-left corner: 3 2 ^ H0(z) z ?1 H^ np?p (z) z ?1 H^ p(z) 6 ^ (z ) H^ q?p (z) z ?1 H^ q+p (z) 777 6 H A^(z) = 66 q.. .. .. 7 5 4 . . . H^ np?q (z) H^ np?q?p (z) H^ np?q+p (z) (assuming q > p).
10
Note that the matrix A^ is no longer Toeplitz. The frequency-domain interpretation of A^ can be obtained by Lemma 1 as follows: v^(z ) = A^(z )^x(z ), where x^ and v^ are the n-fold and (np=q )-fold polyphase vectors of x^(z ) and v^(z ) in Figure 3, respectively. In general, n is a multiple of q . If n = q , to get the matrix A^(z ), we need to perform the (pq )-fold polyphase decomposition of H ; such a fact was also observed by Hsiao [10] from a somewhat di erent point of view. Can we recover H^ (z ) from A^(z )? Since blocking is one-to-one, the equivalent question is: Can we recover H from Dq HEp?
Lemma 5 Assume H is LTI. Then the mapping H 7! Dq HEp is one-to-one if p and q are coprime.
The converse of the lemma is also true, i.e., if p and q are not coprime, then H 7! Dq HEp is not one-to-one.
Proof of Lemma 5 It suces to show that Dq HEp = 0 implies H = 0. So assume
Dq HEp = 0. Let h be the impulse response of H , i.e., h = H, where is the unit impulse signal. It follows that for any integer i, Dq HEpU i = 0, where U is the unit time-delay system. This implies, by the de nition of Ep , Dq HU ip = 0, and, by time invariance of H and the de nition of Dq , h(ip + jq) = 0; 8i; j: (10) Since p and q are coprime, there exist integers p and q such that qp + pq = 1. Thus, for any k, take i = kq and j = kp in (10) to get h(k) = 0; and hence H = 0. 2
Based on Lemma 5, we conclude that the matrix A^(z ) in Theorem 1 contains all (np)-fold polyphase components of H if p and q are coprime.
3.2 Blocking Synthesis Subsystems This subsection relates to blocking the synthesis subsystems Si . Let S be the dual of the system A studied in Section 3.1, namely, the system v 7! y shown in Figure 4, with F LTI and the same assumptions about integers p; q; n: n q > p and n is a multiple of q . We wish to block y by Ln . Since S changes the sampling rates by a factor q=p, we de ne the blocked S to be 1 S = LnSL?np=q : Similarly as in Section 3.1, we can show that S is LTI and the transfer matrix is given by ^ q; S^ = KpFJ 11
v-
-
Eq
-
F
Dp
y-
Figure 4: A synthesis subsystem. where F^ is the (np)-fold polyphase matrix of F and the matrices Jq and Kp are de ned by 82 > > > > 66 > 4 > > : h
2
3
3
9
> 1 1 > > 7 6 > = 0 0 777 7 6 7 6 ; Jq = block .. 7 ; ; 6 .. 7 > 4 . 5 > .5 > > 0 q1 ; 0 q1 h i i Kp = block diag 1 0 0 1p ; ; 1 0 0 1p : Jq has (np=q) diagonal blocks and the size is (np) (np=q); Kp has n diagonal blocks and the size is n (np). In this way we arrive at a dual result.
Theorem 2 Let F^ be the (np)-fold polyphase matrix of F : 2 ^ 3 F (z) z ? F^np? (z) z ? F^ (z) 6 F^ (z) F^ (z) z? F^ (z) 777 6 ^ 6 0
F (z) = 6 4
1
.. .
F^np?1 (z)
1
1
1
1
0
.. .
F^np?2 (z)
1
.. 7: 5 . ^ F0 (z) 2
Then the n (np=q ) transfer matrix for S is obtained by retaining every p-th row and every q-th column of F^ starting from the upper-left corner: 2 ^ 3 F0(z) z ?1 F^np?q (z) z ?1 F^q (z) 6 ^ 7 ?1 ^ ?1 ^ ^S (z ) = 666 Fp.(z ) z Fnp.?q+p (z ) z Fp.+q (z ) 777 .. .. .. 4 5 ^ ^ ^ Fnp?p (z) Fnp?p?q (z) Fq?p (z) (assuming q > p).
Again, in general S^ is no longer Toeplitz but it contains all (np)-fold polyphase components of F if p and q are coprime (Lemma 5). Similar frequency-domain interpretation of S^ can be made by Lemma 1. 12
3.3 Blocking the Overall System Now we shall put together the results in the preceding two subsections to get the n-fold blocking of the overall system T . First, block the analysis and synthesis subsystems Ai and Si in (9) and de ne
Ai = Lnp =q Ai L?n 1 ; Si = LnSi L?np1 =q ; i = 0; 1; ; m ? 1: i
i
i
i
By Theorems 1 and 2, Ai and Si are LTI and their transfer matrices can be computed, respectively, from the (npi )-fold polyphase decompositions of Hi and Fi : A^i is (npi =qi ) n and S^i is n (npi=qi ). Assuming the two matrices A^i and S^i are computed, we proceed to compute T = Ln TL?n 1 . By (8),
T = = =
mX ?1 i=0 mX ?1 i=0 mX ?1 i=0
Ln SiAi L?n 1 (Ln SiL?np1 =q )(Lnp =q Ai L?n 1 ) i
i
i
i
Si Ai:
Hence the transfer matrix for T is
T^(z) =
mX ?1 i=0
S^i (z)A^i (z):
This corresponds to the blocked system in Figure 5, where all systems involved are LTI and signals are blocked, e.g., vi is the (npi =qi )-fold blocking of vi in Figure 1. De ne the analysis matrix and the synthesis matrix to be, respectively, 2 ^ 3 A0(z) 6 ^ 7 h i A^(z) = 664 A1(z) 775 ; S^(z) = S^0(z) S^1 (z) S^m?1 (z) : A^m?1 (z) Note that the analysis (synthesis) matrix consists of polyphase components of the analysis (synthesis) bank only. Thus we conclude
T^(z) = S^(z)A^(z): 13
x
? r
rr
-
A0
-
A1
v0
-
S0
v1
-
S1
vm?1
- Am?
1
- ?j ?r
rr - ?j y -
- Sm?
1
Figure 5: The blocked LTI system. Both the analysis and synthesis matrices, A^ and S^, are n n; for example, since A^i is (npi =qi ) n, the number of columns of A^ is n and the number of rows is
n
mX ?1 p i i=0 qi
= n;
by the assumption in (2). Therefore, the blocked T has a transfer matrix which can be expressed as the product of the synthesis and analysis matrices. In the frequency domain, this means 2 6 4
y^0(z) .. .
y^n?1 (z)
3 2 7 = S^(z )A^(z )6 5 4
x^0 (z) .. .
x^n?1 (z)
3 7; 5
where x^j and y^j are the n-fold polyphase components of the input x and output y in Figure 1 by Lemma 1. Now let us summarize the procedure to compute the n n matrix T^ in Figure 5:
Step 1 For i = 0; 1; ; m ? 1, compute the (npi)-fold polyphase decompositions of Hi and Fi :
H^ i (z) =
npX i ?1 z?j H^ ij (z npi ); j =0
F^i (z) =
npX i ?1 z?j F^ij (znpi ): j =0
Step 2 For i = 0; 1; ; m ? 1, form the (npi)-fold polyphase matrices, H^ i and F^i, associated with Hi and Fi , respectively (Lemma 3); then de ne A^i to be the (npi =qi ) n 14
matrix obtained by retaining every qi -th row and every pi-th column of H^ i , starting from the upper-left corner (Theorem 1), 2 3 H^ i0(z) z ?1 H^ i;np ?p (z) z ?1 H^ ip (z) 6 H^ iq (z) H^ i;q ?p (z) z ?1 H^ i;q +p (z) 777 6 ^ 6 Ai (z) = 6 (11) .. .. .. 7; 4 5 . . . H^ i;np ?q (z) H^ i;np ?q ?p (z) H^ i;np ?q +p (z) i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
i
and similarly, de ne S^i to be the n (npi =qi ) matrix obtained by retaining every pi -th row and every qi -th column of F^i , starting from the upper-left corner (Theorem 2), 2 6 ^Si (z ) = 666 4
F^i0 (z) F^ip (z)
z ?1 F^i;np ?q (z) z ?1 F^iq (z) ? z 1 F^i;np ?q +p (z) z ?1 F^i;p +q (z) i
.. .
i
F^i;np ?p (z) i
.. .
i
i
i
i
i
F^i;np ?p ?q (z) i
i
.. .
i
i
i
i
F^i;q ?p (z) i
3 7 7 7 7: 5
(12)
i
Step 3 Form the two n n analysis and synthesis matrices A^(z) = S^(z) = and then
h h
i0 A^00(z) A^1 0 (z) A^m?1 0(z) ; i S^0(z) S^1 (z) S^m?1 (z)
T^(z) = S^(z)A^(z):
The transfer matrix T^(z ) serves as a frequency-domain model for the general multirate system in Figure 1. The analysis and synthesis matrices A^(z ) and S^(z ) can be regarded as generalization of the polyphase representation matrices for uniform lter banks [29, 26]. Since we assumed that pi and qi are coprime, there exist one-to-one correspondences between A and fHig and between S and fFi g (Lemma 5). In [14] nonuniform lter banks are studied by converting them into uniform lter banks via block diagram manipulations. The frequency-domain model derived here amounts to accomplish the same but uses only matrix manipulations. One can perhaps argue that the method in this paper is more systematical and concise. However, the time-domain approach in [20] is quite di erent. Finally, we remark that a di erent approach for frequency-domain modeling of multirate systems is based on alias-component matrices; this approach was developed for certain multirate systems involving fractional sampling factors in [25, 24]; see also [23]. 15
3.4 Implementation via Blocking We shall conclude Section 3 by looking at an alternative structure, based on the n n analysis and synthesis matrices, to implement the multirate system T in Figure 1. Recall that the cascaded LTI system S A equals the blocked T , namely, Ln TL?n 1 ; hence T = L?n 1 S ALn . It turns out that Ln and L?n 1 can be realized by the n-fold decimator and expander together with the time delay U and advance U .
Lemma 6 The following identities hold: 2 6 6 Ln = 66 4
Dn Dn U .. .
Dn (U )n?1
3 7 7 7 7; 5
i
h
L?n 1 = En UEn U n?1 En :
(13)
The proof is straightforward from the de nitions of the quantities involved. Based on this lemma, we have h
T = En UEn
2 6 i 6 n ? 1 U En S A66 4
Dn DnU .. .
Dn(U )n?1
3 7 7 7 7: 5
This corresponds to the system shown in Figure 6. Note that in this structure, the z transforms of the signals x0 ; x1; ; xn?1 are the n-fold polyphase components of the input x in Figure 1 and similarly for y0 ; y1; ; yn?1 and the output y; the LTI systems A and S operate at a rate which is slower than the normal rate by a factor of n. Such a structure can be used to implement any periodic system based on its blocked, LTI system.
4 Analysis In this section we deal with some analysis issues associated with the multirate system in Figure 1; in particular, we give a solution to the analysis problem stated in Section 1.
4.1 Conditions for Alias Free and Perfect Reconstruction The overall system T in Figure 1 in general su ers aliasing, magnitude distortion, and phase distortion. The system is alias-free if T is LTI and has perfect reconstruction if, in addition, T is free from magnitude and phase distortions; in this case, T equals some pure time-delay 16
x
?
- #n
x0-
y0-
-
"n
z
? r
- #n
x1-
rr
rr r
A
?
y1-
-
z?1
"n
xn?-1
- 6j 6r
rr
S
6
z
- #n
-j y6
yn?1
-
z?1
"n
6
Figure 6: The implementation structure via blocking. system. The conditions for alias free and perfect reconstruction can be stated in terms of the analysis and synthesis matrices, A and S , derived in the preceding section.
Theorem 3 Let T be the nonuniform lter-bank system in Figure 1. (a) T is alias-free i there exists some functions G^ j (z ), j = 0; 1; ; n ? 1, such that 2 ^ 3 G z? G^n? z ? G^ 6 ^ G^ z? G^ 777 ^ ^ 66 G 0
S A = 6 .. 4 . G^ n?1
1
1
1
1
.. .
1
0
G^ n?2
1
.. 7 : . 5 ^ G0 2
(14)
Moreover, if this condition is satis ed, the transfer function for T is
T^(z) =
nX ?1 z?j G^ j (z n ): j =0
(b) T achieves perfect reconstruction i there exists some integers k and r, k 0, 0 r n ? 1, such that " # ?1 Ir 0 z ? k S^ A^ = z ; (15) I 0 n?r
where Ir and In?r are the r r and (n ? r) (n ? r) identity matrices, respectively. Moreover, if this condition is satis ed, T^(z ) = z ?(kn+r) .
17
Proof Part (a) follows directly from Lemma 3. For part (b), recall that T has perfect
reconstruction i T = G, where G is some pure time-delay system: G^ (z ) = z ?d for some d 0. This condition is equivalent to Ln TL?n 1 = LnGL?n 1 , or in the frequency domain, S^(z)A^(z) = G^ (z). To get G^(z), the n-fold polyphase matrix for G, write d = kn + r for some unique integers k and r with k 0 and 0 r n ? 1; then the polyphase components of G are given by ( ?k j = r; ^ Gj (z) = 0z j 6= r: Hence G^ (z ) equals to the right-hand side of (15) by Lemma 3 and the proof is complete. 2 The conditions given in this theorem have been obtained in other forms in [20, 14].
4.2
H1
Analysis
Now we answer the analysis question raised in Section 1 regarding computation of the performance measure J . This is done by evaluating the H1 norm of some transfer matrix. To review, the Hardy space H1 is the Banach space of complex-valued matrix functions of z that are analytic and bounded outside the unit disk, that is, for jz j > 1. Thus H1 is the space of transfer matrices of stable, causal, LTI discrete-time systems. The norm of a matrix G^(z) in H1 is its supremal maximum singular value on the unit circle: h i kG^k1 = sup max G^ ej! : !
If G^ (z ) is the transfer matrix of a stable, causal, LTI discrete-time system G, then the `2 induced norm of G equals the H1 norm of G^ (z ), that is, kGk = kG^ k1 . Back to our problem of computing J . Given the nonuniform system T in Figure 1 and the ideal delay system Tid with T^id (z ) = z ?d for some d 0, the problem is to compute the `2-induced norm in (3), i.e, J = kTid ? T k. This cannot be evaluated directly since T is not LTI. However, blocking Tid ? T preserves the induced norm [Lemma 2 (c)]. So
J = kTid ? T k:
(16)
Now the systems Tid and T are both LTI and so the induced norm in (16) equals the H1 norm of the transfer matrix, J = kT^id ? S^ A^k1; where T^id is derived in Theorem 3 (b): "
0
T^id(z) = z ?k I n?r
#
z ?1 Ir ; d = kn + r; 0 r n ? 1: 0 18
(17)
In summary, given a nonuniform system, computing the reconstruction performance index J is reduced to evaluating the H1 norm of some matrix function involving the analysis and synthesis matrices. The latter problem is standard [7, 8, 4] and MATLAB software exists [1].
5 Design This section studies the optimal design problem stated in Section 1.
5.1
H1
-Optimal Design
For the optimal design problem associated with Figure 1 we assume the analysis lters Hi have been designed based on band-splitting and coding requirements and now we are to design the synthesis lters Fi to minimize J in (3), where T^id (z ) = z ?d and the delay integer d is assumed to be given. Since blocking preserves the induced norm (as in Section 4.2), design of Fi relates to design of the synthesis matrix S ; the equivalent design problem based on the frequency-domain model is now: Given T^id and A, design S^ to minimize J = kT^id ? S^ A^k1 , the H1 norm of the n n transfer matrix T^id ? S^ A^. The matrix T^id depends on the delay integer d and is given in (17); the analysis matrix A^ depends on Hi and is constructed by the procedure in Section 3.3. If S^ can be designed, the synthesis lters Fi are recovered according to the same procedure. The equivalent problem is an H1 optimization problem:
Jopt = inf^ kT^id ? S^ A^k1 : S
(18)
For xed analysis systems, the optimal performance Jopt in general is a function of the delay integer d; so we may write Jopt (d). Computing this function allows us to study the tradeo between reconstruction performance and reconstruction delay. It is a fact proven in [4] for uniform lter banks that lim J (d) = 0: d!1 opt In systems with the same decimation integer, namely, with all equal qi in Figure 1, there is no structural dependency in S^, see Section 6 for a discussion. In this case, we optimize over the class of S^ in (18) which are just stable, causal, and IIR; this is a standard H1 -optimal model-matching problem, which can be solved by standard techniques in control [7, 8, 3] and there exists MATLAB software [1]. For illustration, consider the example below. 19
5.2 Design Example Shown in Figure 7 is a two-channel nonuniform system, in which q0 = q1 = 3, p0 = 2, and
x-
"2 - H
0
- #3
v0-
"3 - F
-
H1
- #3
v1-
"3 - F
- #2
0
? y-j
1
Figure 7: The two-channel example.
p1 = 1. The frequency range [0; ) is split into the low-frequency band [0; 2=3) for channel 0 and the high-frequency band [2=3; ) for channel 1. To get the frequency-domain model, we note that n = 3 and then follow the procedure in Section 3.3. For the analysis matrix, bring in the 6-fold and 3-fold polyphase decompositions of H0 and H1 , respectively,
H^ 0(z) = to get
5 X z?j H^ 0j (z6); j =0
2 A^(z) = 64
H^ 1(z) =
2 X z?j H^ 1j (z3); j =0
3
H^ 00(z) H^ 03(z) H^ 10(z)
z ?1 H^ 04(z) z ?1 H^ 02(z) H^ 01(z) z ?1 H^ 05(z) 75 : z ?1 H^ 12(z) z ?1 H^ 11(z) Similarly, the polyphase decompositions of F0 and F1 , X X F^0(z) = z ?j F^0j (z6); F^1 (z) = z ?j F^0j (z3);
yield the synthesis matrix
5
2
j =0
j =0
2 S^(z) = 64
F^00(z) F^02(z) F^04(z)
z ?1 F^03(z) z ?1 F^05(z) F^01(z)
(19)
3
F^10(z) F^11(z) 75 : F^12(z)
(20)
For an ideal delay system T^id (z ) = z ?d , represent the blocked system Tid in the form of (17) with n = 3. Then the associated H1 optimization problem is Jopt(d) = inf kT^id ? S^ A^k1 : S^
20
If d, H0, and H1 are given, so are T^id and A^; and then the optimization problem can be solved by using the MATLAB function hinfsyn in the -Analysis and Synthesis Toolbox [1]. In this case, the discrete-time problem needs to be converted into an equivalent continuous-time problem via bilinear transformation, see, e.g., [3], and all computations involved are based on state-space data. In this way, the optimal S computed and therefore the optimal F0 and F1 obtained from (20) and (19) are in general IIR. For illustration, we compute the function Jopt (d) for two sets of FIR linear-phase analysis lters, H0 lowpass with cuto frequency =3 (this agrees with the fact that channel 0 processes frequencies up to 2=3 because the 2-fold expander is prior to H0 ) and H1 highpass with cuto frequency 2=3, all obtained by the MATLAB function r1 which uses the Hamming window: Case 1: order of H0 = 19 and order of H1 = 20; Case 2: order of H0 = 23 and order of H1 = 24. The curves for Jopt (d) are shown in Figure 8. The tradeo between reconstruction performance 0
10
−1
10
−2
10
−3
10
−4
10
0
5
10
15
20
25
30
Figure 8: Jopt versus d for Case 1 (solid) and Case 2 (dash). and delay is obvious: For small Jopt , one needs to tolerate relatively large d. It can be proven as in [4] that lim d!1 Jopt (d) = 0 under a mild regularity condition on the analysis matrix, which means that for reasonable analysis lters, one can get arbitrarily close to perfect reconstruction (Jopt = 0) by tolerating suciently large time delay. 21
For a closer look at the optimal synthesis lters computed, we consider the analysis lters in Case 1. The magnitude Bode plots for H0 and H1 are given in Figure 9. Select d = 20. 20
0
−20
−40
−60
−80
−100
−120 0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
Figure 9: Case 1: jH^ 0j (solid) and jH^ 1j (dash) in dB versus !=2 . From the solid curve in Figure 8 we nd Jopt = 0:8066%. The corresponding optimal synthesis matrix S^(z ) is IIR with order 40, from which we compute from (20) and (19) the IIR synthesis lters with their impulse responses shown in Figure 10. The IIR impulse responses, despite their di erent e ective lengths, look very symmetric, indicating that the designed F0 and F1 yield almost linear phase. This suggests to approximate them simply by truncation. Now we truncate the impulse responses, f0 (k) and f1 (k), of F0 and F1 to get their FIR approximations: f0FIR (k) = f0(k); 5 k 56; f1FIR (k) = f1 (k); 0 k 23: With the approximated lters F0FIR and F1FIR , there is a price to be paid in reconstruction performance. Using the analysis result in Section 4.2, the performance achieved by F0FIR and F1FIR is computed to be JFIR = 1:132% (compare with Jopt = 0:8066%). This means that for all input x, the energy gain from x to the reconstruction error e, kek22 =kxk22, is always less than or equal to (1:132%)2 = 0:000128. Figure 11 gives the Bode plots for F0FIR (lowpass) and F1FIR (highpass). With the approximated FIR synthesis lters, the reconstruction errors, e(k) := x(k ? 20) ? y (k), are simulated for x(k) = 1(k), the unit step signal, and x(k) = cos(0:5k) + cos(2:5k), the sum of a low-frequency sinusoid (! = 0:5) and a high-frequency sinusoid (! = 2:5), see Figure 12 for the error signals. 22
3 2 1 0 −1 0
10
20
30
40
50
60
10
20
30
40
50
60
1.5 1 0.5 0 −0.5 −1 0
Figure 10: Impulse responses f0 (k) and f1 (k) versus k for the optimal F0 (upper) and F1 (lower).
20 10 0 −10 −20 −30 −40 −50 −60 −70 −80 0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
Figure 11: jF^0FIRj (solid) and jF^1FIR j (dash) in dB versus !=2 . 23
0.02
0.015
0.01
0.005
0
−0.005
−0.01
−0.015
−0.02 0
10
20
30
40
50
60
70
80
90
100
Figure 12: e(k) versus k for the unit step input (solid) and the sinusoid input (dot). Finally, we remark that better reconstruction performance is possible by using larger time delays, longer synthesis lters, or even other analysis lters.
6 Conclusion The method used in the preceding section for IIR H1 optimal lter design extends directly to the class of nonuniform lter systems with equal decimation integers, namely, the class of systems in Figure 1 with all qi equal. However, if qi are not all the same, there is a diculty. Consider the three-channel case with p0 ; p1 ; p2 = 1 ; 1 ; 1 : q0 q1 q2 2 4 4 In this case n = 4. To write down the synthesis matrix, we introduce the 4-fold polyphase decompositions for Fi , 3 X F^i (z) = z?j F^ij (z 4 ); i = 0; 1; 2; j =0
then the procedure in Section 3.3 yields 2 ^ F00 z?1 F^02 6 ?1 ^ ^ S^ = 664 FF^01 z F^F03 02 00 ^ ^ F03 F01 24
F^10 F^11 F^12 F^13
F^20 3 F^21 777 : F^22 5 F^23
Note that the rst two columns in the synthesis matrix S^ are related to each other: Each is obtained from the four polyphase components of F0 . Such a phenomenon is called structural dependency, which poses a constraint in design. The design approach in Section 5 cannot handle this constraint easily because there the H1 -optimal synthesis matrices are computed from explicit formulas based on solving two Riccati equations [6]. It can be shown that structural dependency occurs whenever qi are not all equal, under the assumption in (2). We remark that this phenomenon was also observed in [14]. One way to handle structural dependency in H1 design is to use numerical optimization. In [22], convex optimization was applied to design H1 -optimal FIR lters in QMF banks. This approach is promising for treating structural dependency in nonuniform lter bank design.
Acknowledgement:
This paper was revised and expanded based on an earlier manuscript submitted to IEEE Trans. on Signal Processing in 1994. The author would like to thank the Associate Editor, Professor Roberto Bamberger, and an anonymous reviewer for encouraging this revision, and to the reviewers for helpful comments.
References [1] G. J. Balas, J. C. Doyle, K. Glover, A. Packard, and R. Smith, -Analysis and Synthesis Toolbox User's Guide, MuSyn Inc. and the MathWorks, Inc., 1991. [2] T. Chen, \A lifting approach to analysis and design of nonuniform multirate lter banks," Proc. Canadian Conf. on Electrical and Computer Engineering, Montreal, Sept. 1995, pp. 874-877. [3] T. Chen and B. A. Francis, Optimal Sampled-Data Control Systems, Springer, London, 1995. [4] T. Chen and B. A. Francis, \Design of multirate lter banks by H1 optimization," IEEE Trans. Signal Processing, vol. 43, pp. 2822-2830, 1995. [5] R. E. Crochiere and L. R. Rabiner, Multirate Digital Signal Processing, Prentice-Hall, Englewood Cli s, NJ, 1983. [6] J. C. Doyle, K. Glover, P. P. Khargonekar, and B. A. Francis, \State-space solutions to standard H2 and H1 problems," IEEE Trans. Automatic Control, vol. 34, no. 8, pp. 831-847, 1989. [7] B. A. Francis, A Course in H1 Control Theory, Springer-Verlag, Lecture Notes in Control and Information Sciences, vol. 88, Berlin, 1987. [8] M. Green and D. J. N. Limebeer, Linear Robust Control, Prentice-Hall, Englewood Cli s, NJ, 1995. 25
[9] P.-Q. Hoang and P. P. Vaidyanathan, \Non-uniform multirate lter banks: theory and design," Proc. IEEE Int. Symp. Circuits and Syst., Portland, OR, 1989, pp. 371-374. [10] C.-C. Hsiao, \Polyphase lter matrix for rational sampling rate conversions," Proc. IEEE Int. Conf. Accoust., Speech, Signal Processing, Dallas, April 1987, pp. 2173-2176. [11] M. R. K. Khansari and A. Leon-Garcia, \Subband decomposition of signals with generalized sampling," IEEE Trans. Signal Processing, vol. 41, no. 12, pp. 3365-3376, 1993. [12] P. P. Khargonekar, K. Poolla, and A. Tannenbaum, \Robust control of linear timeinvariant plants using periodic compensation," IEEE Trans. Automatic Control, vol. 30, pp. 1088-1096, 1985. [13] R. D. Koilpillai and P. P. Vaidyanathan, \A spectral factorization approach to pseudoQMF design," IEEE Trans. Signal Processing, vol. 41, no. 1, pp. 82-92, 1993. [14] J. Kovacevic and M. Vetterli, \Perfect reconstruction lter banks with rational sampling factors," IEEE Trans. Signal Processing, vol. 41, no. 6, pp. 2047-2064, 1993. [15] Y. C. Lim, R. H. Yang, and S.-N. Koh, \The design of weighted minimax quadrature mirror lters," IEEE Trans. Signal Processing, vol. 41, no. 5, pp. 1780-1789, 1993. [16] B. Liu and L. T. Bruton, \The design of nonuniform-band maximally decimated lter banks," Proc. IEEE Int. Symp. Circuits and Syst., Chicago, IL, 1993, pp. 375-378. [17] R. A. Meyer and C. S. Burrus, \A uni ed analysis of multirate and periodically timevarying digital lters," IEEE Trans. Circuits and Systems, vol. 22, pp. 162-168, 1975. [18] S. Mirabbasi, B. A. Francis, and T. Chen, \Controlling distortions in maximally decimated lter banks," IEEE Trans. on Circuits and Systems II: Analog and Digital Signal Processing, to appear. [19] T. Miyawaki and C. W. Barnes, \Multirate recursive digital lters { a general approach and block structures," IEEE Trans. Acoustics, Speech, Signal Processing, vol. 31, pp. 1148-1154, 1983. [20] K. Nayebi, T. P. Barnwell III, and M. J. T. Smith, \Nonuniform lter banks: a reconstruction and design theory," IEEE Trans. Signal Processing, vol. 41, no. 3, pp. 1114-1127, 1993. [21] T. Q. Nguyen, \Near-perfect-reconstruction pseudo-QMF banks," IEEE Trans. Signal Processing, vol. 42, no. 1, pp. 65-76, 1994. [22] S. Ratzla , An Optimal Model Matching Approach to Multirate FIR Filter Bank Design, M.Sc. Thesis, Dept. of Electrical and Computer Engineering, University of Calgary, August 1995. [23] R. L. Reng, \Polyphase and modulation descriptions of multirate systems { a systematic approach," Proc. Int. Conf. DSP, 1995, pp. 212-217. 26
[24] R. G. Shenoy, \Analysis of multirate components and application to multirate lter design," Proc. IEEE Int. Conf. Accoust., Speech, Signal Processing, 1994, vol. 3, pp. 121-124. [25] R. G. Shenoy, D. Burnside, and T. W. Parks, \Linear periodic systems and multirate lter design," IEEE Trans. Signal Processing, vol. 42, pp. 2242-2256, 1994. [26] P. P. Vaidyanathan, \Theory and design of M -channel maximally decimated Quadrature mirror lters with arbitrary M , having the perfect-reconstruction property," IEEE Trans. Acoustics, Speech, Signal Processing, vol. 35, no. 4, pp. 476-492, 1987. [27] P. P. Vaidyanathan, Multirate Systems and Filter Banks, Prentice-Hall, Englewood Cli s, NJ, 1993. [28] P. P. Vaidyanathan and S. K. Mitra, \Polyphase networks, block digital ltering, LPTV systems, and alias-free QMF banks: a uni ed approach based on pseudocirculants," IEEE Trans. Acoustics, Speech, Signal Processing, vol. 36, no. 3, pp. 381-391, 1988. [29] M. Vetterli, \A theory of multirate lter banks," IEEE Trans. Acoustics, Speech, Signal Processing, vol. 35, no. 3, pp. 356-372, 1987.
27
