- Email: [email protected]
Recursive calculation of expansion coefficients for generalized orthonormal rational bases
14th World Congress ofIFAC
RECURSIVE CALCULATION OF EXPANSION COEFFICIENTS FO...
H-3a-II-4
Copyright © 1999 lFAC
14th Triennial World Congress» Beijing, P.R. China
RECURSIVE CALCULATION OF EXPANSION COEFFICIENTS FOR GENERALIZED ORTHONORMAL RATIONAL BASES Per Bodin
S3-A1ttomatic Control The Royal Institute of Technology 100 44 Stockholm;r Sweden Email: [email protected]
Abstract: A recursive scheme is derived for calculation of the expansion co€fficients of a known system transfer function into a generalized orthonormal rational basis. The inherent structure of the basis functions leads to a computationally efficient method, which nee.ds the solution of a low-dimensional Sylvester equation a.t each step. The method is demonstrated \vith a brief example. Possible applications are also discussed. Copyright © 1999lFAC Keywords: All pass filters, Filter banks, Function approximation~ Least-squares approximation, R.ecursive algorithms, State-space realization, Transfer functions
1.
I~TRODuCTIOj\T
The possibility of using orthonormal rational basis functions in system identification has become a field of considerable interest over the last few years. The properties of the classical first order Lagucrre and second order Kautz constructions were examined for system identification in (\~lahlberg;, 1991; '~lahlberg, 1994). These bases were generalized in (Heubergcr et al., 1995) to a const.ruction based on the repetition of an inner function of any order.
In (l'iinlless and Gustafsson, 1997; Bodin ct al., 1996), even more general constructions were given. (Bodin et al., 1996) gives a historical overvie\v of similar constructions concluding that some of these results were known by Takenaka and by IVlalmquist already in 1925; see (Takenaka, 1925; Walsh~ 1960). It also gives a convenient way, based on the results of (Roberts and Mullis, 1987), to construct basis functions. Basically, such a basis can be constructed from an array of serially connected all-pass filters where each filter has a balanced state space realization. The ba..~is functions are then given by the transfer functions from
the input to the array to the states of each of the filters. These filters are allu\\red to he different and even of different order. al.~ 1998; Bodin et al.~ 1997; Bodin, 1997), the observation that a different basis for the same span is obtained if the filters in the array are rearranged leads to a scheme for selecting the best basis with respect to some cost function. From ideas similar to those for the best \vavelet packet selection of (Coifman and Wickerhauscr, 1992), the best orthonormal rational £1_ or entropy-basis can be chosen. In many cases, this leads to compact expansions concentrating energy into only a few large coefficients. The observation that the number of different basis functions is smaller than the nUlnber of possible basis constructions gives algorithms evaluated in a number of steps of a complexity no larger than the number of basis functions. This situation is sinlilar to that of wavelet packets methods. The algorithm however, becomes tot.ally different. If for example, there are m different all-pass filters, the best of m~ different bases is selected in only O(m2-m-l) steps.
In (Bodin et
4141
Copyright 1999 IFAC
ISBN: 0 08 043248 4
RECURSIVE CALCULATION OF EXPANSION COEFFICIENTS FO...
14th World Congress of IFAC
It is straightforward to sho,,,", that the filter Hk is orthogonal all-pass if and only if
The best basis selection needs the coefficients of all possible basis functions as input. In (Bodin et al., 1998; Badin and \Vahlberg, 1998; Bodin, 1997) a scheme is derived that, given the expansion coefficients in one basis, calculates the coefficients for all the other possible basis functions in the same complexity as the best basis selection. This is achieved using local orthogonal matrix transformations.The kno~7n, or given coefficients can for example be obtained as estimates from system identification; see for example (Bodin and Wahlberg, 1994).
F
[~: ~:]
:=
(4)
=
is an orthogonal matrix, i.e. FF T (Roberts and Mullis, 1987).
I; see
Define a vector valued basis function as the transfer function between the input u to the array and the state outputs Xk of each of the filters H k ~ k = 1) .. ~ ~ m. lvIore specifically, define a basis function as
In this paper, a scheme is given for the calculation of expansion coefficients for a known sy-stem transfer functioIl. The rnethod is recursive and perInits coefficient calculation only needing the solution of a v x nk Sylvcster equation in each step, where v and nk are the respective McMillan degrees of the kno\vn system and the basis all-pass filter at the current step. The total number of steps is equal to the total number of all-pass filters.
and
In Section 2, the construction of generalized orthonormal rational bases is briefly described. The
The 1-l 2 -inner product betvi,"een tions f and h is defined as
(5) where qlk(Z)
coefficient calculation scheme is derived in Section 3 and a brief example is given in Section 4. Possible applications are discussed in Section 5 and finally in Section 6, some conclusions are made.
= (zI - Ak)-l B
(6)
k
is the McMillan degree of H k •
nk
t\VQ
scalar func-
J 7r
1
{f, h) := 2 7f
f(eiw)h(e-iw)dw,
(7)
and similarly, the outer product (f, h)o between two column vector functions f and h can be
vlritten (f,h)o :=
The perhaps simplest ~"ay to construct a generalized orthonormal rational basis was revealed by the results of (Roberts and Mullis, 1987). Consider an array of serially connected all-pass filters as shown in Figure 1, where each of the filters has a balanced state space realization. In (Roberts and Mullis, 1987), such a filter is referred to as an orthonormal all-pass filter.
1L(t)~~
--~ Xl
~
... --~
(zI - A
x~
k )
~1 Bk
+ Dk,
(8)
f(eiw)hT(e-iw)ckv.
\vhere (A ul , B:rn) is the state space realization of 'l'(z). Then, it is easy to sho\v that ('It, lft) 0 === I if
ArnA~ + Bn1:B~ == I.
(10)
• Consider two orthogonal all-pass filters Hi and H2 "\vith
1t1ore specifically, let the filter H k (z) have the state space realization (A k , B k , C k , Dk) so that
= Ck
21f
• Let
Fig. 1. An array of serially connected all-pass filters with input balanced realizations~
Hk(z)
~
Orthonormality of the family {llrk}k=l is then fulfilled for t"vo reasons:
~
X2
J 'r.
2. PR,ELIMINARIES
F1
(1)
v'lhere
:=
[~~ ~~]
and
F 2 :=
[~~ ~~]
and respective states Xl and X2. T'he serial connection HI followed by H 2 with the combined state
(2)
and (11)
(3)
4142
Copyright 1999 IFAC
ISBN: 0 08 043248 4
RECURSIVE CALCULATION OF EXPANSION COEFFICIENTS FO...
14th World Congress of IFAC
The expansion coefficient gk is given by
is also orthogonal and all-pass. This follows from the fact that the serial connection has
F;:=
[
Al
0
Bl]
B 2 C 1 A 2 B 2D 1 D2C I C 2 D 2D 1
vlhich is orthogonal if Fl and F anal.
2
g'[ ;:::
are orthogg:=
[:J
(18)
can in the same way be written
gT = (G, 'ls)o
It \vas sho",'n already in (Takenaka~ 1925) that if the array is expanded so that rn -7 cx>, the
(19)
where 'IJ(z) is defined by (9) _ The coefficients g can thus be expressed as
functions {Wk}~l are complete in 1{2 (omitting exactly proper functions) if and only if
gT == cP,
nk
L L 1-
IPk~ll ==
00,
(13)
(20)
v,rith
k=11=1
1r
p;= 2~J(e'.wI-a)-lbB~(e-iWI~A~)-ldwl
\vhere {Pk,I}~~l arc the poles of Hk. This quite mild condition simply states that the poles of the all-pass filters must not (if at all) converge to the unit circle too fast as m -+ 00.
(21) where (.Am., Bm-, C-rn, Drn ) is the state space realization for the combined state
1\lany special cases can now easily be obtained from this way of constructing rational bases. If all the filters {Hk}k=l are equal, the construction of (Heuberger et al., 1995) is obtained. _~ll-pass filters of second order gives the Kautz functions while first order filters lead to the Laguerre basis; see (,~rahlberg, 1991; \Vahlberg, 1994)~
(22) and where P solves the discrete-time Sylvester equation
N O\V, given a scalar strictly proper transfer function G (z ) with all of its poles inside the unlt circle. Let n ::=: L::~=1 nk and define the nth order approximation G n of G as Gn(Z) :=
(17)
a.ccording to the notation established in the previous section. All of the coefficients
(12)
These two facts imply orthonormality of the functions {Wk lr::=l"
0::
(G, '1' k>O
L
gZ\Pk(Z).
aPA~
::=
(23)
P.
Ho\vever, for large n this eqnation may be difficult, or inefficient to solve~ In order to \vork around this problcm~ the inherent structure of the basis funetions can be ernployed. T'his leads to a scheme derived be}ol,v, only needing the solution of a v x nk-dimensional Sylv€ster equation in each of the k == 1, ... ,1T~ recursive steps.
(14)
k=l
It is well known that the mean squared error GnIJ~ is minimized if
+ bB~
IIC -
Given the m all-pass filters {Hk } k=l and their balanced state space realizations (Ak,Bk,Ck,Dk) for k ~ 1, ... , rn. The state space realization (Ak,Bk,Ck,D k ) of the kth combined state
gk,l ] [ (G, tPk,l >] gk := : = ,: = {G, '1'd~. (15) [ gk,'nk (G, V'k,nk) In (Bodin and Wahlberg, 1994), it was demonstrated how the coefficients {gk}r:::l could be estimated from noisy input/output data from G using system identification and thresholding. The following section gives a recursive scheme for calculating {gk}~=l in the case of G being a known transfer function.
(24)
(25)
Ch
3. COEFFICIENT CALCUL..t\.TION
for k Assume that the strictly proper and stable transfer function G(z) is known and has 1v1cMillan degree v. Suppose that it has a minima.l state space realization Ca, b, c) such that
G(z) = c (zI - a)-l b.
Let
~==
==
Pk
[DkC k -
1
Ck]
Dk;:=
DkDk-l,
2, ... ,m with the initial conditions
Al
:=
CJ
:==
AI, Cl,
Rt
:=
D1
:= D 1 -
Ht,
(26)
solve
(27)
(16)
4143
Copyright 1999 IFAC
ISBN: 0 08 043248 4
RECURSIVE CALCULATION OF EXPANSION COEFFICIENTS FO...
Then, according to (20) and
[~:]
g[:=
T
14th World Congress of IFAC
which has one pole in 0.7 and two poles in 0.25 ± O.15i. Select the minimal state space realization
(21)~
=
1.2 -0.2175 2 0 [ o 0.25 c == [0 0.5 0.2] .
(28)
cPko
a==
Introduce the partitioning
corresponding to the partit.ioning in (25). The Sylvester equation (27) then becomes T
-
-T
T aPkAk+aP k _ 1 Ck_1B k
T
-
+bDk_1B k
(i) : 1 pole in 0.5 (3 filters) (ii) : 3 poles in 0.2 and 0.3 ± O.2i (1 filter) (iii) : 2 poles in 0.4 ± O.3i (2 filters).
(30) (31)
~Pk.
In this V\ray~ there is a total of 10 states in the basis. The ordering of the filters in the array was {(i), (ii), (iii), (i), (iii), (i)}. The coefficients calculated using Algorithm 1 arc shown in Figure 2 (a). For COInparison, Figure 2 Cb) sho\vs the least squares coefficient estimate from a noisefree simulation of (41). Finally, Figure 2 Cc) shows
Define
Vk :==
aPkCf + bDk.
(32)
From (25) and (29), V k follows the recursion V k = Vk-1Dk + aPkCf VI = aP 1 + bD1 .
er
0 0
The basis was generated from 6 all-pass filters of three different types
(29)
aP k-l A[-1 + bB[_l = P k-l
00119]
(33) (34)
the magnitude of the difference between these t,vn ca.."ies.
Vlith the definitions above, the following algorithm can no",,~ be formed:
Algorithm 1 ~ (Recursive coefficient calculation).
5. APPLICATIONS
1. Let PI be the solution of
aP1Ai'
+ bBi' = PI
Except from the immediate possibility of analyzing a known~ or identified transfer function in an orthonormal fixed-pole expansion, one other possibility is to use the algoritlun introduced above for modeling additive disturbances in best basis selection. The complete framevlork presented in (Bodin et al., 1997; Bodin and 'VVahlberg, 1998; Bodin et al., 1998) permits best basis selection with respect to a large variety of different cost functions such as the entropy- or el-costs. Onc common property for these is that the best basis should provide as compact an expansion a.s possible. In this ~vay, there will be many small coefficients but only a few large.
(35)
and let
g'f := cP 1 VI := PI
er + bD
(36) 1.
(37)
2. For k = 2, ... ,m, solve the Sylvester equation
and let
g'f :== cP k' As long as k
V
k
(39)
< rn, let also :== Vk~lDk
+aPkC[.
(40)
From the derivation above, if follo\vs that all the coefficients {gk } k~ 1 are calculated by this scheme. One possible "\vay to solve the Sylvestcr equation is suggested in Appendix A. 'The solution is based on Kronecker products. The following section gives an example of coefficient calculation using Algorithm 1.
4. EXAl\-1PLE
Consider the given system G( ) z
+ G.1 + 0.085) (z
z = (Z2 _ O.5z
- 0.7)
(41)
In (Saito and Coifman, 1995), Local DiscriIninant. Basis (LDB) selection was considered for different libraries of bases such as ,vavelet packets or local cosines. In these methods, the best basis is selected virith respect to a modified, or weighted cost function taking into account the coefficient distribution of for example an additive disturbance model. \Veighting is made so that it is expensive to select coordinates "",~here disturbances are kno"\\,"'n to be large, but inexpensive to select a coordinate where almost no disturbances at all are present. In this way an orthonormal basis is selected as ttfar away as possible'~ from the disturbance~ or from some other class of functions. In (Borlin and Villemoes, 1997), this ,vas proposed with a new cost function for speech enhancement. The scheme introduced in this paper permits coefficient calculation for the transfer function of some 4144
Copyright 1999 IFAC
ISBN: 0 08 043248 4
14th World Congress ofIFAC
RECURSIVE CALCULATION OF EXPANSION COEFFICIENTS FO...
known model of an additive disturbance. T'he coefficients can then be calculated for all different filter orderings using the methods of (Bodin and '"Vahlberg, 1998) and thus incorporated into the cost function in order to achieve LDD selection also fOT gener alized orthonormaI rational bases.
Exact ExpansIon Coefficients
1.5
6. CONCLUDING DISCUSSION The recursive scheme presented in this paper gives a convenient way to calculate the expansion of a known system int.o a generalized orthonormal ratioIlal basis. The algorithm ho\vever, needs to be further analyzed.
0.5
o
2
3
4 5 6 7 8 Coefficient Number
9
10
Firstly, the numerical condition of the sohltion to the discrete-time Sylvester equation (38) has not been considered here. Its dependence on the chosen state space realization of the known system (16) has to be clarified. The expansion coefficients themselves are obviously independent of the choice of such a realization. However, the solution of (38) is not~ so guidelines for selecting t.he realization are clearly needed.
Ca) Expansion coefficients generated llsing Algorithm 1.
Estimated Expansion Coefficients
1.5
Secondly, since the coefficients are independent of the realization of the known system, it would be desirable to obtain a scheme for calculating coefficients which is also realization independent. It remains to shoV\,r if this is feasible.
0.5
o
2
3
4
5
6
7
8
9
As mentioned in Section 5, an obvious application of Algorithm 1 lies within Local Discriminant Basis selection for generalized orthonormal rational bases. This is examined in (Bodin, 1999).
10
Coefficient Nu mber (b) Least squares estimate of the expansion coefficients.
X
Bodin~
10-4 Coefficient Error Magnitude
1
2
3
4 5 6 7 8 Coefficient Numbe r
7. REFERENCES
9
10
(c) Magnitude of the coefficient error betVlo'een the t",'o cases.
Fig. 2. Coefficients calculated in the example.
P. (1997). On the Selection of Best Orthonormal Basis in System Identification and Signal Processing. PhD thesis. KTH. Stockholm, S)veden. TRIT..~-REG-9702. Bodin, P. (1999). Selection of local discriminant generalized orthonormal rational basis. In: Proc. 14th IFAG World Con.qress. Beijing, China. To appear. Bodin, P. and B. Wahlberg (1994). Thresholding in high order transfer function estimation. In: Pr'oc. 3.'rrd IEEE Conf. on Decision and Control. ,rol. 4. Orlando, FL. pp. 3400-3405. Bodin, P. and B. vVahlberg (1998). Fast expansion coefficient calculation for best orthonormal rational ba.~is_ In: Proc. MTNS'98. Padova, Italy. To appear. Bodin, P. and L. F. 'lillemoes (1997). Spectral subtraction in the time-frequency domain using wavelet packets. In: Proc. IEEE Workshop on Speech Coding- Pocono Manor, Pennsylvania. pp. 47-48. 4145
Copyright 1999 IFAC
ISBN: 008 0432484
RECURSIVE CALCULATION OF EXPANSION COEFFICIENTS FO...
14th World Congress of IFAC
Bodin, P., L. F. Villemoes and B. ';V"'ahlberg (1997). An algorithm for selection of best orthonormal rational basis. In: Proc. 86th IEEE Conf. on Decision and Control. San Diego) CA. pp. 1277-1282. Bodin, P., L. F. Villemoes and B. Wahlberg (1998). Selection of best orthonormal rational basis. SIAM Journal on Control and Optimization. Accepted. Bodin, P., T. Oliveira e Silva and B. Wahlberg (1996). On the construction of orthonormal basis functions for system identification. In: Proc. 13th IFAC World Congress. \Tol. 1. San Francisco~ CA. pp. 369-374. Coifman, R. R. and ::VI. V. \Vickerhauser (1992). Entropy-ba..~ed algorithms for best basis selection. IEEE Trans. on Information Theory 38(2), 713-718. Heuberger, P. S. C.~ P. M. J. Van den Hof and O. H. Bosgra (1995). A generalized orthonormal basis for linear dynamical systems. IEEE Transactions on Automatic Control 40(3), 451-465. Ninness~ B. M. and F. Gustafsson (1997). A unifying construction of orthonormal bases for system identification. IEEE Transactions on A utomatic Control 42(4), 515-521. Roberts, R~ A. and C. T. Ivlullis (1987). Digital Signal Processing. Addison-Wesley Series in Electrical Engineering: Digital Signal Processing. Addison-\V'"esley Publishing Cornpany. Reading, :f\,lA. Saito~ N. and R. R. Coifman (1995). Local discriminant bases and their applications. Journal of Mathematical Imaging and Vision 5(4),337-358. Takenaka, S. (1925). On the orthogonal functions and a new formula of interpolation. Japanese Journal of Mathematics 11, 129-145. 'v~ahlberg, B~ (1991). System identification using Laguerre models. IEEE Transactions on Automatic Control 36(5), 551-562. W~ahlberg, B. (1994). System identification using Kautz models. IEEE Transactions on Auto~ matic Control 39(6), 1276-1281. \\7alsh, J. L. (1960). Interpolation and Approximation by Rational Functions in the Complex Domain. VaL XX of American Mathematical Society Colloquium Publications~ third ed .. American Ivlathematical Society. Providence, Rhode Island. First edition in 1935.
v..There A l and A 2 are n'l x nl and n2 x n2 matrices respectively and where B I and B2 are nl x 1 and 1~2 x 1 column vectors respectively. Let
where Pk are nl x 1 column vectors. Define and
"vhere Then
I~
denotes the Kronccker tensor product.
-
P :=
Pt] [ :
-
1-
== (I - A)- B.
(A.4)
Pn2
Appendix ..t\.. SOLUTION OF THE SYLVESTER EQUATION Given the equation A1PAr
+ B 1 B;f
:::=
P,
(A.I) 4146
Copyright 1999 IFAC
ISBN: 0 08 043248 4
RECURSIVE CALCULATION OF EXPANSION COEFFICIENTS FO...
H-3a-II-4
Copyright © 1999 lFAC
14th Triennial World Congress» Beijing, P.R. China
RECURSIVE CALCULATION OF EXPANSION COEFFICIENTS FOR GENERALIZED ORTHONORMAL RATIONAL BASES Per Bodin
S3-A1ttomatic Control The Royal Institute of Technology 100 44 Stockholm;r Sweden Email: [email protected]
Abstract: A recursive scheme is derived for calculation of the expansion co€fficients of a known system transfer function into a generalized orthonormal rational basis. The inherent structure of the basis functions leads to a computationally efficient method, which nee.ds the solution of a low-dimensional Sylvester equation a.t each step. The method is demonstrated \vith a brief example. Possible applications are also discussed. Copyright © 1999lFAC Keywords: All pass filters, Filter banks, Function approximation~ Least-squares approximation, R.ecursive algorithms, State-space realization, Transfer functions
1.
I~TRODuCTIOj\T
The possibility of using orthonormal rational basis functions in system identification has become a field of considerable interest over the last few years. The properties of the classical first order Lagucrre and second order Kautz constructions were examined for system identification in (\~lahlberg;, 1991; '~lahlberg, 1994). These bases were generalized in (Heubergcr et al., 1995) to a const.ruction based on the repetition of an inner function of any order.
In (l'iinlless and Gustafsson, 1997; Bodin ct al., 1996), even more general constructions were given. (Bodin et al., 1996) gives a historical overvie\v of similar constructions concluding that some of these results were known by Takenaka and by IVlalmquist already in 1925; see (Takenaka, 1925; Walsh~ 1960). It also gives a convenient way, based on the results of (Roberts and Mullis, 1987), to construct basis functions. Basically, such a basis can be constructed from an array of serially connected all-pass filters where each filter has a balanced state space realization. The ba..~is functions are then given by the transfer functions from
the input to the array to the states of each of the filters. These filters are allu\\red to he different and even of different order. al.~ 1998; Bodin et al.~ 1997; Bodin, 1997), the observation that a different basis for the same span is obtained if the filters in the array are rearranged leads to a scheme for selecting the best basis with respect to some cost function. From ideas similar to those for the best \vavelet packet selection of (Coifman and Wickerhauscr, 1992), the best orthonormal rational £1_ or entropy-basis can be chosen. In many cases, this leads to compact expansions concentrating energy into only a few large coefficients. The observation that the number of different basis functions is smaller than the nUlnber of possible basis constructions gives algorithms evaluated in a number of steps of a complexity no larger than the number of basis functions. This situation is sinlilar to that of wavelet packets methods. The algorithm however, becomes tot.ally different. If for example, there are m different all-pass filters, the best of m~ different bases is selected in only O(m2-m-l) steps.
In (Bodin et
4141
Copyright 1999 IFAC
ISBN: 0 08 043248 4
RECURSIVE CALCULATION OF EXPANSION COEFFICIENTS FO...
14th World Congress of IFAC
It is straightforward to sho,,,", that the filter Hk is orthogonal all-pass if and only if
The best basis selection needs the coefficients of all possible basis functions as input. In (Bodin et al., 1998; Badin and \Vahlberg, 1998; Bodin, 1997) a scheme is derived that, given the expansion coefficients in one basis, calculates the coefficients for all the other possible basis functions in the same complexity as the best basis selection. This is achieved using local orthogonal matrix transformations.The kno~7n, or given coefficients can for example be obtained as estimates from system identification; see for example (Bodin and Wahlberg, 1994).
F
[~: ~:]
:=
(4)
=
is an orthogonal matrix, i.e. FF T (Roberts and Mullis, 1987).
I; see
Define a vector valued basis function as the transfer function between the input u to the array and the state outputs Xk of each of the filters H k ~ k = 1) .. ~ ~ m. lvIore specifically, define a basis function as
In this paper, a scheme is given for the calculation of expansion coefficients for a known sy-stem transfer functioIl. The rnethod is recursive and perInits coefficient calculation only needing the solution of a v x nk Sylvcster equation in each step, where v and nk are the respective McMillan degrees of the kno\vn system and the basis all-pass filter at the current step. The total number of steps is equal to the total number of all-pass filters.
and
In Section 2, the construction of generalized orthonormal rational bases is briefly described. The
The 1-l 2 -inner product betvi,"een tions f and h is defined as
(5) where qlk(Z)
coefficient calculation scheme is derived in Section 3 and a brief example is given in Section 4. Possible applications are discussed in Section 5 and finally in Section 6, some conclusions are made.
= (zI - Ak)-l B
(6)
k
is the McMillan degree of H k •
nk
t\VQ
scalar func-
J 7r
1
{f, h) := 2 7f
f(eiw)h(e-iw)dw,
(7)
and similarly, the outer product (f, h)o between two column vector functions f and h can be
vlritten (f,h)o :=
The perhaps simplest ~"ay to construct a generalized orthonormal rational basis was revealed by the results of (Roberts and Mullis, 1987). Consider an array of serially connected all-pass filters as shown in Figure 1, where each of the filters has a balanced state space realization. In (Roberts and Mullis, 1987), such a filter is referred to as an orthonormal all-pass filter.
1L(t)~~
--~ Xl
~
... --~
(zI - A
x~
k )
~1 Bk
+ Dk,
(8)
f(eiw)hT(e-iw)ckv.
\vhere (A ul , B:rn) is the state space realization of 'l'(z). Then, it is easy to sho\v that ('It, lft) 0 === I if
ArnA~ + Bn1:B~ == I.
(10)
• Consider two orthogonal all-pass filters Hi and H2 "\vith
1t1ore specifically, let the filter H k (z) have the state space realization (A k , B k , C k , Dk) so that
= Ck
21f
• Let
Fig. 1. An array of serially connected all-pass filters with input balanced realizations~
Hk(z)
~
Orthonormality of the family {llrk}k=l is then fulfilled for t"vo reasons:
~
X2
J 'r.
2. PR,ELIMINARIES
F1
(1)
v'lhere
:=
[~~ ~~]
and
F 2 :=
[~~ ~~]
and respective states Xl and X2. T'he serial connection HI followed by H 2 with the combined state
(2)
and (11)
(3)
4142
Copyright 1999 IFAC
ISBN: 0 08 043248 4
RECURSIVE CALCULATION OF EXPANSION COEFFICIENTS FO...
14th World Congress of IFAC
The expansion coefficient gk is given by
is also orthogonal and all-pass. This follows from the fact that the serial connection has
F;:=
[
Al
0
Bl]
B 2 C 1 A 2 B 2D 1 D2C I C 2 D 2D 1
vlhich is orthogonal if Fl and F anal.
2
g'[ ;:::
are orthogg:=
[:J
(18)
can in the same way be written
gT = (G, 'ls)o
It \vas sho",'n already in (Takenaka~ 1925) that if the array is expanded so that rn -7 cx>, the
(19)
where 'IJ(z) is defined by (9) _ The coefficients g can thus be expressed as
functions {Wk}~l are complete in 1{2 (omitting exactly proper functions) if and only if
gT == cP,
nk
L L 1-
IPk~ll ==
00,
(13)
(20)
v,rith
k=11=1
1r
p;= 2~J(e'.wI-a)-lbB~(e-iWI~A~)-ldwl
\vhere {Pk,I}~~l arc the poles of Hk. This quite mild condition simply states that the poles of the all-pass filters must not (if at all) converge to the unit circle too fast as m -+ 00.
(21) where (.Am., Bm-, C-rn, Drn ) is the state space realization for the combined state
1\lany special cases can now easily be obtained from this way of constructing rational bases. If all the filters {Hk}k=l are equal, the construction of (Heuberger et al., 1995) is obtained. _~ll-pass filters of second order gives the Kautz functions while first order filters lead to the Laguerre basis; see (,~rahlberg, 1991; \Vahlberg, 1994)~
(22) and where P solves the discrete-time Sylvester equation
N O\V, given a scalar strictly proper transfer function G (z ) with all of its poles inside the unlt circle. Let n ::=: L::~=1 nk and define the nth order approximation G n of G as Gn(Z) :=
(17)
a.ccording to the notation established in the previous section. All of the coefficients
(12)
These two facts imply orthonormality of the functions {Wk lr::=l"
0::
(G, '1' k>O
L
gZ\Pk(Z).
aPA~
::=
(23)
P.
Ho\vever, for large n this eqnation may be difficult, or inefficient to solve~ In order to \vork around this problcm~ the inherent structure of the basis funetions can be ernployed. T'his leads to a scheme derived be}ol,v, only needing the solution of a v x nk-dimensional Sylv€ster equation in each of the k == 1, ... ,1T~ recursive steps.
(14)
k=l
It is well known that the mean squared error GnIJ~ is minimized if
+ bB~
IIC -
Given the m all-pass filters {Hk } k=l and their balanced state space realizations (Ak,Bk,Ck,Dk) for k ~ 1, ... , rn. The state space realization (Ak,Bk,Ck,D k ) of the kth combined state
gk,l ] [ (G, tPk,l >] gk := : = ,: = {G, '1'd~. (15) [ gk,'nk (G, V'k,nk) In (Bodin and Wahlberg, 1994), it was demonstrated how the coefficients {gk}r:::l could be estimated from noisy input/output data from G using system identification and thresholding. The following section gives a recursive scheme for calculating {gk}~=l in the case of G being a known transfer function.
(24)
(25)
Ch
3. COEFFICIENT CALCUL..t\.TION
for k Assume that the strictly proper and stable transfer function G(z) is known and has 1v1cMillan degree v. Suppose that it has a minima.l state space realization Ca, b, c) such that
G(z) = c (zI - a)-l b.
Let
~==
==
Pk
[DkC k -
1
Ck]
Dk;:=
DkDk-l,
2, ... ,m with the initial conditions
Al
:=
CJ
:==
AI, Cl,
Rt
:=
D1
:= D 1 -
Ht,
(26)
solve
(27)
(16)
4143
Copyright 1999 IFAC
ISBN: 0 08 043248 4
RECURSIVE CALCULATION OF EXPANSION COEFFICIENTS FO...
Then, according to (20) and
[~:]
g[:=
T
14th World Congress of IFAC
which has one pole in 0.7 and two poles in 0.25 ± O.15i. Select the minimal state space realization
(21)~
=
1.2 -0.2175 2 0 [ o 0.25 c == [0 0.5 0.2] .
(28)
cPko
a==
Introduce the partitioning
corresponding to the partit.ioning in (25). The Sylvester equation (27) then becomes T
-
-T
T aPkAk+aP k _ 1 Ck_1B k
T
-
+bDk_1B k
(i) : 1 pole in 0.5 (3 filters) (ii) : 3 poles in 0.2 and 0.3 ± O.2i (1 filter) (iii) : 2 poles in 0.4 ± O.3i (2 filters).
(30) (31)
~Pk.
In this V\ray~ there is a total of 10 states in the basis. The ordering of the filters in the array was {(i), (ii), (iii), (i), (iii), (i)}. The coefficients calculated using Algorithm 1 arc shown in Figure 2 (a). For COInparison, Figure 2 Cb) sho\vs the least squares coefficient estimate from a noisefree simulation of (41). Finally, Figure 2 Cc) shows
Define
Vk :==
aPkCf + bDk.
(32)
From (25) and (29), V k follows the recursion V k = Vk-1Dk + aPkCf VI = aP 1 + bD1 .
er
0 0
The basis was generated from 6 all-pass filters of three different types
(29)
aP k-l A[-1 + bB[_l = P k-l
00119]
(33) (34)
the magnitude of the difference between these t,vn ca.."ies.
Vlith the definitions above, the following algorithm can no",,~ be formed:
Algorithm 1 ~ (Recursive coefficient calculation).
5. APPLICATIONS
1. Let PI be the solution of
aP1Ai'
+ bBi' = PI
Except from the immediate possibility of analyzing a known~ or identified transfer function in an orthonormal fixed-pole expansion, one other possibility is to use the algoritlun introduced above for modeling additive disturbances in best basis selection. The complete framevlork presented in (Bodin et al., 1997; Bodin and 'VVahlberg, 1998; Bodin et al., 1998) permits best basis selection with respect to a large variety of different cost functions such as the entropy- or el-costs. Onc common property for these is that the best basis should provide as compact an expansion a.s possible. In this ~vay, there will be many small coefficients but only a few large.
(35)
and let
g'f := cP 1 VI := PI
er + bD
(36) 1.
(37)
2. For k = 2, ... ,m, solve the Sylvester equation
and let
g'f :== cP k' As long as k
V
k
(39)
< rn, let also :== Vk~lDk
+aPkC[.
(40)
From the derivation above, if follo\vs that all the coefficients {gk } k~ 1 are calculated by this scheme. One possible "\vay to solve the Sylvestcr equation is suggested in Appendix A. 'The solution is based on Kronecker products. The following section gives an example of coefficient calculation using Algorithm 1.
4. EXAl\-1PLE
Consider the given system G( ) z
+ G.1 + 0.085) (z
z = (Z2 _ O.5z
- 0.7)
(41)
In (Saito and Coifman, 1995), Local DiscriIninant. Basis (LDB) selection was considered for different libraries of bases such as ,vavelet packets or local cosines. In these methods, the best basis is selected virith respect to a modified, or weighted cost function taking into account the coefficient distribution of for example an additive disturbance model. \Veighting is made so that it is expensive to select coordinates "",~here disturbances are kno"\\,"'n to be large, but inexpensive to select a coordinate where almost no disturbances at all are present. In this way an orthonormal basis is selected as ttfar away as possible'~ from the disturbance~ or from some other class of functions. In (Borlin and Villemoes, 1997), this ,vas proposed with a new cost function for speech enhancement. The scheme introduced in this paper permits coefficient calculation for the transfer function of some 4144
Copyright 1999 IFAC
ISBN: 0 08 043248 4
14th World Congress ofIFAC
RECURSIVE CALCULATION OF EXPANSION COEFFICIENTS FO...
known model of an additive disturbance. T'he coefficients can then be calculated for all different filter orderings using the methods of (Bodin and '"Vahlberg, 1998) and thus incorporated into the cost function in order to achieve LDD selection also fOT gener alized orthonormaI rational bases.
Exact ExpansIon Coefficients
1.5
6. CONCLUDING DISCUSSION The recursive scheme presented in this paper gives a convenient way to calculate the expansion of a known system int.o a generalized orthonormal ratioIlal basis. The algorithm ho\vever, needs to be further analyzed.
0.5
o
2
3
4 5 6 7 8 Coefficient Number
9
10
Firstly, the numerical condition of the sohltion to the discrete-time Sylvester equation (38) has not been considered here. Its dependence on the chosen state space realization of the known system (16) has to be clarified. The expansion coefficients themselves are obviously independent of the choice of such a realization. However, the solution of (38) is not~ so guidelines for selecting t.he realization are clearly needed.
Ca) Expansion coefficients generated llsing Algorithm 1.
Estimated Expansion Coefficients
1.5
Secondly, since the coefficients are independent of the realization of the known system, it would be desirable to obtain a scheme for calculating coefficients which is also realization independent. It remains to shoV\,r if this is feasible.
0.5
o
2
3
4
5
6
7
8
9
As mentioned in Section 5, an obvious application of Algorithm 1 lies within Local Discriminant Basis selection for generalized orthonormal rational bases. This is examined in (Bodin, 1999).
10
Coefficient Nu mber (b) Least squares estimate of the expansion coefficients.
X
Bodin~
10-4 Coefficient Error Magnitude
1
2
3
4 5 6 7 8 Coefficient Numbe r
7. REFERENCES
9
10
(c) Magnitude of the coefficient error betVlo'een the t",'o cases.
Fig. 2. Coefficients calculated in the example.
P. (1997). On the Selection of Best Orthonormal Basis in System Identification and Signal Processing. PhD thesis. KTH. Stockholm, S)veden. TRIT..~-REG-9702. Bodin, P. (1999). Selection of local discriminant generalized orthonormal rational basis. In: Proc. 14th IFAG World Con.qress. Beijing, China. To appear. Bodin, P. and B. Wahlberg (1994). Thresholding in high order transfer function estimation. In: Pr'oc. 3.'rrd IEEE Conf. on Decision and Control. ,rol. 4. Orlando, FL. pp. 3400-3405. Bodin, P. and B. vVahlberg (1998). Fast expansion coefficient calculation for best orthonormal rational ba.~is_ In: Proc. MTNS'98. Padova, Italy. To appear. Bodin, P. and L. F. 'lillemoes (1997). Spectral subtraction in the time-frequency domain using wavelet packets. In: Proc. IEEE Workshop on Speech Coding- Pocono Manor, Pennsylvania. pp. 47-48. 4145
Copyright 1999 IFAC
ISBN: 008 0432484
RECURSIVE CALCULATION OF EXPANSION COEFFICIENTS FO...
14th World Congress of IFAC
Bodin, P., L. F. Villemoes and B. ';V"'ahlberg (1997). An algorithm for selection of best orthonormal rational basis. In: Proc. 86th IEEE Conf. on Decision and Control. San Diego) CA. pp. 1277-1282. Bodin, P., L. F. Villemoes and B. Wahlberg (1998). Selection of best orthonormal rational basis. SIAM Journal on Control and Optimization. Accepted. Bodin, P., T. Oliveira e Silva and B. Wahlberg (1996). On the construction of orthonormal basis functions for system identification. In: Proc. 13th IFAC World Congress. \Tol. 1. San Francisco~ CA. pp. 369-374. Coifman, R. R. and ::VI. V. \Vickerhauser (1992). Entropy-ba..~ed algorithms for best basis selection. IEEE Trans. on Information Theory 38(2), 713-718. Heuberger, P. S. C.~ P. M. J. Van den Hof and O. H. Bosgra (1995). A generalized orthonormal basis for linear dynamical systems. IEEE Transactions on Automatic Control 40(3), 451-465. Ninness~ B. M. and F. Gustafsson (1997). A unifying construction of orthonormal bases for system identification. IEEE Transactions on A utomatic Control 42(4), 515-521. Roberts, R~ A. and C. T. Ivlullis (1987). Digital Signal Processing. Addison-Wesley Series in Electrical Engineering: Digital Signal Processing. Addison-\V'"esley Publishing Cornpany. Reading, :f\,lA. Saito~ N. and R. R. Coifman (1995). Local discriminant bases and their applications. Journal of Mathematical Imaging and Vision 5(4),337-358. Takenaka, S. (1925). On the orthogonal functions and a new formula of interpolation. Japanese Journal of Mathematics 11, 129-145. 'v~ahlberg, B~ (1991). System identification using Laguerre models. IEEE Transactions on Automatic Control 36(5), 551-562. W~ahlberg, B. (1994). System identification using Kautz models. IEEE Transactions on Auto~ matic Control 39(6), 1276-1281. \\7alsh, J. L. (1960). Interpolation and Approximation by Rational Functions in the Complex Domain. VaL XX of American Mathematical Society Colloquium Publications~ third ed .. American Ivlathematical Society. Providence, Rhode Island. First edition in 1935.
v..There A l and A 2 are n'l x nl and n2 x n2 matrices respectively and where B I and B2 are nl x 1 and 1~2 x 1 column vectors respectively. Let
where Pk are nl x 1 column vectors. Define and
"vhere Then
I~
denotes the Kronccker tensor product.
-
P :=
Pt] [ :
-
1-
== (I - A)- B.
(A.4)
Pn2
Appendix ..t\.. SOLUTION OF THE SYLVESTER EQUATION Given the equation A1PAr
+ B 1 B;f
:::=
P,
(A.I) 4146
Copyright 1999 IFAC
ISBN: 0 08 043248 4