- Email: [email protected]
Frequency-domain model reduction approach to design IIR digital filters using orthonormal bases
Int. J. Electron. Commun. (AEÜ) 60 (2006) 413 – 420 www.elsevier.de/aeue
Frequency-domain model reduction approach to design IIR digital filters using orthonormal bases Rabah W. Aldhaheri∗ Department of Electrical and Computer Engineering, King Abdulaziz University, P.O. Box 80204, Jeddah 21589, Saudi Arabia Received 8 August 2004; received in revised form 20 May 2005
Abstract A new and numerically efficient technique for the design of reduced-order IIR digital filters is presented. This technique is based on computing the cross-gramian matrix from the state–space representation of the full-order filter in an arbitrary frequency range of operation. The modified cross-gramian is derived in the frequency domain and computed by solving only one Lyapunov equation. This technique does not require a minimal system to start with, and the reduced order can be obtained without calculating the conventional balancing transformation. Instead, orthonormal bases are computed to find the left and right eigenspaces associated with the large eigenvalues of the modified cross-gramian matrix. These eigenspaces are used to convert the high-order to a reduced-order models. The phase linearity is also discussed and it is shown that the proposed method can be used to transform a linear-phase FIR filter to a reduced-order IIR filter and the phase linearity is preserved over the interested frequency band, which is commonly be the passband. A comparison between the conventional time-domain methods such as balanced truncation (BT) and optimal Hankel-norm approximation (OHA); and the proposed frequency-domain model reduction technique is presented through several filters design examples to illustrate the effectiveness of the proposed technique. 䉷 2005 Elsevier GmbH. All rights reserved. Keywords: IIR filters; Frequency-domain cross-gramian; Model reduction; Filter approximation
1. Introduction In many engineering applications such as signal processing, estimation/filtering and control, it is necessary to approximate the high-order complex system by low-order simple model. The resulting low-order model is less expensive and more efficient computationally and hence it is used as a replacement to the original high-order model. Over the last few decades, a great deal of attention has been paid to model reduction techniques and many model reduction algorithms
∗ Tel.: +966 2 695 2195; fax: +966 2 695 2686.
E-mail address: [email protected]. 1434-8411/$ - see front matter 䉷 2005 Elsevier GmbH. All rights reserved. doi:10.1016/j.aeue.2005.05.016
have been proposed. Earlier methods were mainly based on retaining the dominant poles in the reduced model as in singular perturbation approximation [1], or on matching some moments of the original system as in Padé approximation techniques [2]. An interesting approach for model reduction is the one in which each component of the state vector has equal information about controllability and observability, i.e., the input-to-state coupling and the state-to-output coupling are weighted equally. This approach is known in the literature as the balanced state–space realization [3–7]. The reducedorder is determined by truncating the system states to those associated with the most significant singular values and the remaining states are rejected. The drawback of balanced
414
Rabah W. Aldhaheri / Int. J. Electron. Commun. (AEÜ) 60 (2006) 413 – 420
realization-based algorithms is that it requires the original system to be completely controllable and observable and consequently, these algorithms are ill-conditioned when the original system is close to being non- minimal. In [8,9], an equivalent to the balanced reduced model is proposed without computing the balancing transformation. Moreover, it was shown in [9] that the balanced reduced model can be found even when the original model is non-minimal. Other interesting approaches for model reduction are the ones based on Hankel-norm optimal approximation [10,11], and impulse response gramians [12]. In digital filter design, the aforementioned model reduction algorithms are used to approximate the high-order filters (often FIR) to reduced-order IIR filters [6,7,11–14]. The major reason for starting with FIR filters is because the requirements of linear phase can always be met in the design. Then the FIR is approximated by IIR filters and the linearity of the phase can be preserved, at least, in the passband region. In these techniques, the maximum magnitude error between the full-order and the reduced order is bounded from above by twice the sum of the singular values associated with the eliminated states. Hence, there is still need for improving the error of the reduced filters. In [13], the singular perturbation technique is shown to improve the error in the low-frequency range, and in [14], the finite and optimal search methods are presented to reduce the error in the passband. In [15], the frequencydomain-balanced structure approach is developed and it was shown that the balanced realization obtained by this approach shares many similar properties with its time-domain counterpart. In this paper, we extend the balanced structure approach and introduce the concept of frequency-domain cross-gramian. Moreover, a modified Lyapunov equation is derived to solve the frequency-domain cross-gramian matrix over any specified frequency range, the band where we would like to reduce the approximation error over it. The eigenspace of this matrix is used to convert the large-order filters to low-order filters. Unlike the method in [15], the proposed method does not solve two Lyapunov equations and avoids computing the balanced transformation. This paper is organized as follows: In Section 2, we review some aspects of balanced discrete systems and give the necessary background to formulate the problem addressed in this paper. In Section 3, the concept of the frequencydomain cross-gramian matrix and how it can be computed is presented. In Section 4, we show how to compute orthogonal bases which span the right eigenspaces associated with the large eigenvalues of the frequency-domain cross-gramian matrix using numerically stable method. The Sylvester equation is used to compute the left eigenspaces. The proposed algorithm is also summarized in this section. Examples of digital filters design are given in Section 5 to illustrate the significance of the proposed method and conclusions are given in Section 6.
2. Preliminaries and problem formulation Consider an nth-order, linear, time-invariant and asymptotically stable discrete-time system described by X(k + 1) = AX(k) + Bu(k), y(k) = CX(k) + Du(k),
(2.1)
where A ∈ Rn×n , B ∈ Rn×1 and C ∈ R1×n . The controllability and observability gramians are defined as follows: WC = WO =
∞ k=0 ∞
Ak BB T (AT )k , (AT )k C T CAk .
(2.2)
k=0
These gramians satisfy the following Lyapunov equations: AW C AT − WC + BB T = 0,
(2.3a)
AT WO A − WO + C T C = 0.
(2.3b)
Fernando and Nicholson [16] demonstrated that the information in the controllability gramian and the observability gramian are essentially contained in the cross-gramian matrix defined by WCO =
∞
Ak BCAk .
(2.4)
k=0
Equivalently, WCO can be computed by solving the Lyapunov equation AW CO A − WCO + BC = 0.
(2.5)
We refer to the controllability, observability and crossgramian defined by (2.2)–(2.5) as the time-domain gramians in contrast to the frequency-domain gramians, which will be defined later in following section. The controllability and the observability gramians are used in many previous works [3–5] to get a balanced system realization {Ab , Bb , Cb } provided that the original system is minimal. This was achieved by finding a transformation matrix, T such that 1 0 . (2.6) T −1 WC T −T = T T WO T = = 0 2 The entries of the diagonal of 1 and 2 are assumed to be in a descending order and 1 = diag(1 , 2 , . . . , r ) and 2 = diag(r+1 , r+2 , . . . , n ).
(2.7)
The balanced realization system is determined by A11 A12 B1 , Bb = T −1 B = Ab = T −1 AT = A21 A22 B2 and Cb = CT = [ C1 C2 ] , (2.8)
Rabah W. Aldhaheri / Int. J. Electron. Commun. (AEÜ) 60 (2006) 413 – 420
where A11 ∈ Rr×r , B1 ∈ Rr×1 and C1 ∈ R1×r . The subsystem {A11 , B1 , C1 , D} is considered as a good reducedorder approximation of the original system if r ?r+1 . This rth-order subsystem is called a direct truncation approximation of the balanced system or shortly balanced truncation (BT) approximation. In [8], an equivalent reduced-order model to the subsystem {A11 , B1 , C1 } is proposed without calculating the balanced transformation matrix, but with calculating a transformation matrix that decomposes the original system into strongly controllable and observable coupled subsystem and weakly coupled subsystem. The strongly coupled is retained and the weakly coupled is discarded. Moreover in [13], the reduced-order singularly perturbed system is calculated to obtain a better error approximation in the low-frequency range. These algorithms are numerically more robust compared with the algorithms based on balanced transformation. Yet, they do not give good approximation behavior over the desired frequency range. This motivates us to look for a new technique to achieve a better approximation in any frequency range that we are interested in. In the next section, we propose a definition for the frequency-domain cross-gramian and we will show how it can be used for getting a reduced-order model which behaves closely enough to the original model in the desired frequency range. The model that we consider in this paper is a model of digital filters of different types.
3. Frequency-domain cross-gramian In this section, the frequency-domain cross-gramian matrix is defined and the modified Lyapunov equation to compute this matrix is derived. This definition is an extension to the frequency-domain controllability and observability gramians given in [15]. Moreover, the reduced-order technique is different and it is not based on balanced realization as we will see later in the following section. Lemma 3.1. For the discrete-time system described by (2.1), the cross-gramian matrix is given by 1 WCO = (ej I − A)−1 BC(e−j I − A)−1 d, (3.1) 2 − where I is an n × n identity matrix. Proof. The proof is a direct application of Parseval’s theorem to Eq. (2.4).
415
is equal to the time-domain cross-gramian represented by Eq. (2.4). Theorem 3.1. The frequency-domain cross-gramian matrix, WF , defined by (3.2) satisfies the Lyapunov equation AW F A − WF = −KBC − BCK ∗ , where
(3.3)
∗
denotes the complex conjugate and 2 1 K= (I − e−j A)−1 d − I, 2 1 4
(3.4)
where = 2 − 1 radians/second. Proof. Adding and subtracting the terms ej WCO A, e−j AW CO and WCO to Eq. (2.5) to obtain − (ej I − A)WCO (e−j I − A) + WCO (e−j I − A)ej + (ej I − A)WCO e−j = BC. (3.5) Multiply Eq. (3.5) from the left by (ej I − A)−1 , and from the right by (e−j I − A)−1 . Then, take the integration 1/2 12 d for both sides to obtain WF = −
WCO + K1 WCO + WCO K1∗ , 2
(3.6)
where
2 1 K1 = (I − e−j A)−1 d. 2 1 Substituting (3.6) into the LHS of (3.3), yields ∗ A − WCO + K1 WCO + WCO K1 A 2 − − WCO + K1 WCO + WCO K1∗ . 2
(3.7)
(3.8)
It is easy to show that AK 1 = K1 A and AK ∗1 = K1∗ A. Therefore, Eq. (3.8) can be rewritten as −
[AW CO A − WCO ] + K1 [AW CO A − WCO ] 2 + [AW CO A − WCO ]K1∗ = − K1 − I BC − BC K1∗ − I 4 4 ∗ = −KBC − BCK = RHS of (3.3).
It can be shown that if 1 = − and 2 = , then the Lyapunov equation of the frequency-domain cross-gramian defined by (3.3) is equivalent to its time-domain counterpart defined by (2.5).
Definition 3.1. The frequency-domain cross-gramian in the frequency range (1 , 2 ) is defined as 2 1 WF = (ej I − A)−1 BC(e−j I − A)−1 d. (3.2) 2 1
4. Model reduction algorithm using the corresponding eigenspaces of the frequency-domain cross-gramian
If the whole frequency range, i.e., (−, ) is selected in the above integration, the frequency-domain cross-gramian
In this section, we consider a transformation that decomposes the cross-gramian matrix, WF (or WCO ) into two
416
Rabah W. Aldhaheri / Int. J. Electron. Commun. (AEÜ) 60 (2006) 413 – 420
blocks, one with large eigenvalues and the other with small eigenvalues. This transformation separates the strongly coupled controllable and observable states from the weakly coupled controllable and observable states based on the large and small eigenvalues of the cross-gramian matrix. This transformation block diagonalizes the matrix WF into two blocks: one with large eigenvalues and the other with small eigenvalues, as follows: S 0 V −1 WF V = S = L , (4.1) 0 SS where SL and SS are r × r and (n − r) × (n − r) matrices, respectively. i = |i (SL )|, i = 1, 2, . . . , r, are the large (in magnitude) eigenvalues of WF , and i = |i (SS )|, i = r + 1, r + 2, . . . , n are the small (in magnitude) eigenvalues of WF . Now, let the matrices V and V −1 be partitioned as UL −1 , (4.2) V = [ VL VS ] and V = US where UL [ VL US
I VS ] = r 0
0 In−r
.
VS ] = [ C 1
C2 ]
and
D = D.
(4.4)
The reduced-order model is defined as {A11 , B 1 , C 1 , D} and the transfer function of the reduced model is defined as Hr (z) = D + C 1 [zI r − A11 ]−1 B 1 .
W˜ 11 P − P W˜ 22 + W˜ 12 = 0.
(4.7)
Notice here, the Sylvester equation is computed for two upper triangular forms, W˜ 11 and W˜ 22 as a byproduct of the ordered real Schur form. This will reduce the computational effort required for computing the matrix P . Ir P Step 4: The matrix block diagonalizes W˜ F in 0 In−r
(4.6) and puts it in the form of (4.1). Thus if we partition Q in (4.6) as Q=[ Q1 Q2 ], where Q1 ∈ Rn×r , then VL =Q1 and UL = QT1 − P QT2 . Step 5: The reduced-order model state–space representation is defined as 1. 2. 3. 4.
A11 = UL AV L B 1 = UL B C 1 = CV L D = D.
(4.3)
The columns of VL and the rows of UL span the right and the left eigenspaces associated with (SL ), i.e., the large eigenvalues of WF . Similarly, the columns of VS and the rows of US span the right and the left eigenspaces associated with (SS ), i.e., the small eigenvalues of WF . Applying the transformation, Eq. (4.2), to the original system (2.1), we obtain UL A11 A12 A[ VL VS ] = A = , US A21 A22 UL B1 B= , B= US B2 C = C[ VL
where the eigenvalues of WF appear in descending order of absolute value along the diagonal of W˜ F , W˜ 11 and W˜ 22 are r × r and (n − r) × (n − r) matrices, respectively. Step 3: Compute the matrix P ∈ Rr×(n−r) by solving the Sylvester equation
(4.5)
In the rest of this section, let us summarize the proposed algorithm to obtain a reduced-order model and in the next section we apply this algorithm to design filter-order reduction in an arbitrary frequency range. Step 1: Given the original state–space representation {A, B, C, D} of the full-order digital filter, determine the frequency range, (1 , 2 ) over which it is desired to approximate the full-order filter by a reduced-order one. Then, compute WF either by solving the Lyapunov equation as in Eq. (3.3), or equivalently, by numerical integration of (3.2). Step 2: Compute the ordered real Schur form of WF [8,17]. W˜ 11 W˜ 12 QT WF Q = W˜ F = , (4.6) 0 W˜ 22
Substitute into (4.5) to obtain the transfer function of the reduced-order filter Hr (z).
5. Numerical examples In this section, we present three examples of filter-order reduction using the proposed method and two of the widely used time-domain-based algorithms: the BT and the optimal Hankel-norm approximation (OHA) methods to illustrate the effectiveness of the proposed method in obtaining lower magnitude error in the desired frequency range. Our concern in the first two examples is the magnitude response and how to get a better error performance. In the third example the phase response is discussed and it is shown that the phase linearity can always be preserved provided that the full-order filter is of linear phase. The design procedure consists of two steps: first we design a filter that meets certain specifications using one of the standard procedures. Second, we approximate the full-order filter with a reducedorder one using the algorithm proposed in Section 4 such that the aforementioned specifications are met, at least, in the chosen frequency band. Example 1. Design a Butterworth lowpass that meets the following specifications: 1, 0 0.45, j |H (e )| = 0, 0.54 . The minimum stopband attenuation is 25 dB. It is found that an 12th-order IIR filter meets the above specifications. The magnitude response, which we will
Rabah W. Aldhaheri / Int. J. Electron. Commun. (AEÜ) 60 (2006) 413 – 420 -3
1.4 4
Full order, n = 12 BT method, r = 7 OHA method, r = 7 Proposed method, r = 7
1.2
x 10
3 2 Magnitude response of error
1 Magnitude response
417
0.8
0.6
0.4
1 0 -1 -2 -3
0.2
BT method, r = 7 OHA method, r = 7 Proposed method, r = 7
-4
0
0
0.1
0.2
0.3
0.4 0.5 0.6 Normalized frequency
0.7
0.8
0.9
-5
1
Fig. 1. Magnitude responses of the original and reduced-order IIR filters.
0.2
0.3 0.4 Normalized frequency
0.5
0.6
0.7
1.4
x 10
BT method, r = 7 OHA method, r = 7 Proposed method, r = 7
4.5
0.1
Fig. 3. Magnitude responses error of the two reduced-order filters.
-3
5
0
Full order, n = 6 BT method, r = 4 OHA method, r = 4 Proposed method, r = 4
1.2
4
1 Magnitude response
Error magnitude
3.5 3 2.5 2
0.8
0.6
0.4
1.5 1
0.2
0.5 0
0 0
0.1
0.2
0.3 0.4 Normalized frequency
0.5
0.6
0.7
0
0.1
0.2
0.3
0.4 0.5 0.6 Normalized frequency
0.7
0.8
0.9
1
Fig. 2. Error magnitude of the reduced-order IIR filters.
Fig. 4. Magnitude responses of the original and reduced-order highpass IIR filters.
consider it here as the full order, is depicted in Fig. 1. By applying the algorithm proposed in Section 4, a 7th-order IIR filter is selected as a good approximate to the full-order filter over the frequency band [−0.55, 0.55]. Moreover, a 7th-order IIR filter is computed using the conventional BT and OHA methods. Fig. 1 shows the magnitude responses of the full order (dashed curves), the proposed reduced order (solid curves), and the BT method (dotted curves), and the OHA (dashed-dotted curves) are plotted over the normalized frequency, where the normalized frequency 1 represents
the Nyquist frequency. As we notice, the magnitude response of the proposed method is much closer to the full order than the reduced order obtained by the BT and OHA methods. Figs. 2 and 3 show error magnitude and magnitude response error over the normalized frequency band [0, 0.55]. These figures show clearly that the error over the selected band for the proposed method is much better than that for the BT and OHA methods.
H (z) =
Example 2. Consider the highpass IIR filter represented by the transfer function
4.7297 × 10−2 (1 − 6z−1 + 15z−2 − 20z−3 + 15z−4 − 6z−5 + z−6 ) . 1 − 0.486z−1 + 1.2052z−2 − 0.042145z−3 + 0.36827z−4 + 0.08425z−5 + 0.04481z−6 The magnitude response of this filter is shown in Fig. 4. Again, we applied the BT and OHA methods, and the
418
Rabah W. Aldhaheri / Int. J. Electron. Commun. (AEÜ) 60 (2006) 413 – 420 0.08
10
BT method, r = 4 OHA method, r = 4 Proposed method, r = 4
FIR filter, n = 45 BT method, r = 15 OHA method, r = 15 Proposed method, r = 15
0
0.04
-10
0.02
-20
Magnitude in dB
Magnitude response error
0.06
0 -0.02
-30 -40 -50
-0.04
-60
-0.06
-70
-0.08 0.3
0.4
0.5
0.6 0.7 Normalized frequency
0.8
0.9
1
Fig. 5. Magnitude response error of the three reduced-order filters.
0
0.1
0.2
0.3
0.4 0.5 0.6 Normalized frequency
0.7
0.8
0.9
1
Fig. 6. Magnitude responses of lowpass FIR and reduced-order IIR filters.
proposed method with selected frequency range [0.35, ] to obtain a 4th-order IIR filters represented, respectively, by HBT (z) =
10−2 (4.7297 − 27.255z−1 + 64.4025z−2 − 71.306z−3 + 38.592z−4 ) , 1 − 0.28218z−1 + 0.80547z−2 + 0.10036z−3 + 0.14751z−4
HOHA (z) =
10−2 (4.7297 − 26.684z−1 + 63.728z−2 − 69.916z−3 + 38.805z−4 ) 1 − 0.26310z−1 + 0.84277z−2 + 0.10496z−3 + 0.19029z−4
and Hr (z) =
10−2 (4.7297 − 27.307z−1 + 72.732z−2 − 81.597z−3 + 56.266z−4 ) . 1 − 0.53423z−1 + 0.83175z−2 − 0.032746z−3 + 0.10221z−4
The magnitude responses of the lower-order filters obtained by the BT, OHA and the proposed methods are shown in Fig. 4. To compare between the three approximations, the magnitude response error of the BT method and OHA, which is based on time-domain model reduction, and the proposed method, which is based on frequency-domain model reduction, are plotted in Fig. 5. It is shown that the proposed method gives better error performance than that of the BT and OHA methods in the passband [0.4,1]. It is worth mentioning here that the order reduction ratio from FIR to IIR filters is much more than the reduction ratio when we approximate the IIR by another IIR filters. This will be addressed in the following example. Example 3. Design a linear-phase lowpass FIR filter with the following specifications: j
|H (e )| =
1, 0,
0 0.3, 0.35 ,
where the maximum passband attenuation is 1 dB and the minimum stopband attenuation is less than 25 dB. Using the Remez algorithm [18], it is found that an FIR of order 45 meets the above requirements. The magnitude response of this filter is shown in Fig. 6. Now, if we apply the proposed algorithm with a selected frequency band [−0.36, 0.36], an IIR filter of order 15 is a good approximation which tracks the FIR responses fairly closely over the selected band. The magnitude responses of the lowerorder filters obtained by the BT, OHA and the proposed method are shown in Fig. 6. A comparison between the three approximations is depicted in Fig. 7. As it is clearly shown in this figure, the proposed method gives a better error performance over the selected passband [0 − 0.36] compared with that of the BT and OHA methods, respectively. Fig. 8 shows the corresponding group delay for both the FIR filter and the three reduced-order approximations. Notice that the phase obtained by the proposed algorithm is linear over the band [0 − 0.36]. This confirms that the proposed algorithm maintains the phase linearity as it is the case with the time-domain-based algorithms.
Rabah W. Aldhaheri / Int. J. Electron. Commun. (AEÜ) 60 (2006) 413 – 420
phase linearity is preserved over the selected band provided that the original filter has a linear phase.
0.025 BT method, r = 15 OHA method, r = 15 Proposed method, r = 15
0.02
419
Magnitude response error
0.015 0.01
References
0.005 0 -0.005 -0.01 -0.015 -0.02 -0.025
0
0.05
0.1
0.15 0.2 0.25 Normalized frequency
0.3
0.35
0.4
Fig. 7. Magnitude responses error of the three approximated IIR filters.
60 FIR filter, n = 45 BT method, r = 15 OHA method, r = 15 Proposed method, r = 15
40
Group delay
20
0
-20
-40
-60
0
0.1
0.2
0.3
0.4 0.5 0.6 Normalized frequency
0.7
0.8
0.9
1
Fig. 8. Group delay responses for the FIR and the three approximated IIR filters.
6. Conclusions In this paper, the frequency-domain cross-gramian, WF is derived and computed by solving a Lyapunov equation (3.3). The left and the right eigenspaces associated with the large eigenvalues of WF are used to transform the full-order to a reduced-order models without computing the balancing transformation which tends to be ill-conditioned, especially when the original model is non-minimal or when it has very nearly uncontrollable or unobservable modes. This technique is applied to several digital filters designs and it is found that the error performance of the filters design by this technique outperforms the time-domain-based algorithms such as the balanced truncation method and the optimal Hankel-norm approximation over the selected frequency range. Moreover, using the proposed technique, the
[1] Kokotovic P, Khalil H, O’Reilly J. Singular perturbation methods in control: analysis and design. New York: Academic Press; 1986. [2] Shamash Y. Model reduction using the Routh stability criterion and Pade approximations. Int J Control 1975;21: 475–84. [3] Moore BC. Principal component analysis in linear systems: controllability, observability and model reduction. IEEE Trans Automat Control 1981;AC-26:17–32. [4] Pernebo L, Silverman LM. Model reduction via balanced state–space representation. IEEE Trans Automat Control 1982;AC-27:382–7. [5] Fernando KV, Nicholson H. On the structure of balanced and other principal representations of SISO systems. IEEE Trans Automat Control 1983;AC-28:228–31. [6] Beliczynski B, Kale I, Cain GD. Approximation of FIR by IIR digital filters: an algorithm based on balanced model reduction. IEEE Trans Signal Process 1992;40:532–42. [7] Kale I, Gryka J, Cain GD, Beliczynski B. FIR filter order reduction: balanced model truncation and Hankel-norm optimal approximation. IEE Proc Vision Image Signal Process 1994;141:168–74. [8] Aldhaheri RW. Model order reduction via real Schur form decomposition. Int J Control 1991;53:709–16. [9] Aldhaheri RW, Bamani A. Balanced realization of nonminimal MIMO discrete time systems. Int J Systems Sci 1994;25:173–82. [10] Glover K. All optimal Hankel norm approximation of linear multivariable systems and their L∞ error bounds. Int J Control 1984;39:1115–93. [11] Chen B-S, Peng S-C, Chiou B-W. IIR filter design via optimal Hankel norm approximation. IEE Proc G Circuits Devices Systems 1992;139:586–90. [12] Sreeram V, Agathoklis P. Design of linear-phase IIR filters via impulse-response gramians. IEEE Trans Signal Process 1992;40:389–94. [13] Aldhaheri RW. Design of linear-phase IIR digital filters using singular perturbational model reduction. IEE Proc Vision Image Signal Process 2000;147:409–14. [14] Li L, Xie L, Yan W-Y, Soh YC. Design of low-order linearphase IIR filters via orthogonal projection. IEEE Trans Signal Process 1999;47:448–57. [15] Wang D, Zilouchian A. Model reduction of discrete linear systems via frequency-domain balanced structure. IEEE Trans Circuits Systems 2000;47:830–7. [16] Fernando KV, Nicholson H. On the cross-gramians for symmetric MIMO systems. IEEE Trans Circuits Systems 1985;CAS32:487–9. [17] Laub AJ. A Schur method for solving algebraic Riccati equations. IEEE Trans Automat Control 1979;AC-24:913–21. [18] Proakis JG, Manolakis DG. Digital signal processing: principles, algorithms, and applications. 3rd ed., Englewood Cliffs, NJ: Prentice-Hall; 1996.
420
Rabah W. Aldhaheri / Int. J. Electron. Commun. (AEÜ) 60 (2006) 413 – 420
Rabah W. Aldhaheri was born in Jeddah, Saudi Arabia in 1953. He received his B.S. degree from Riyadh (currently, King Saud) University in 1976, the M.S. degree from Ohio University, in 1981, and the Ph.D. degree from Michigan State University in 1988, all in Electrical Engineering. From 1976 to 1981, he was with the P.T.T. Ministry, Jeddah, Saudi Arabia as an Electrical Engineer and then as a
Director of the Microwave and Multiplex Stations. In 1981, he joined King Abdulaziz University as a Lecturer in the Department of Electrical and Computer Engineering. He currently is an Associate Professor in the same department. Dr. Aldhaheri has held visiting research scholar positions with Michigan State University during 1994–1995, and Queensland University of Technology (QUT) in Brisbane, Australia in 2000. His research interests include digital signal processing with application to filter design, speaker recognition, face recognition and iris recognition; and digital communications. He is an Associate Editor of the Journal of Applied Sciences.
Frequency-domain model reduction approach to design IIR digital filters using orthonormal bases Rabah W. Aldhaheri∗ Department of Electrical and Computer Engineering, King Abdulaziz University, P.O. Box 80204, Jeddah 21589, Saudi Arabia Received 8 August 2004; received in revised form 20 May 2005
Abstract A new and numerically efficient technique for the design of reduced-order IIR digital filters is presented. This technique is based on computing the cross-gramian matrix from the state–space representation of the full-order filter in an arbitrary frequency range of operation. The modified cross-gramian is derived in the frequency domain and computed by solving only one Lyapunov equation. This technique does not require a minimal system to start with, and the reduced order can be obtained without calculating the conventional balancing transformation. Instead, orthonormal bases are computed to find the left and right eigenspaces associated with the large eigenvalues of the modified cross-gramian matrix. These eigenspaces are used to convert the high-order to a reduced-order models. The phase linearity is also discussed and it is shown that the proposed method can be used to transform a linear-phase FIR filter to a reduced-order IIR filter and the phase linearity is preserved over the interested frequency band, which is commonly be the passband. A comparison between the conventional time-domain methods such as balanced truncation (BT) and optimal Hankel-norm approximation (OHA); and the proposed frequency-domain model reduction technique is presented through several filters design examples to illustrate the effectiveness of the proposed technique. 䉷 2005 Elsevier GmbH. All rights reserved. Keywords: IIR filters; Frequency-domain cross-gramian; Model reduction; Filter approximation
1. Introduction In many engineering applications such as signal processing, estimation/filtering and control, it is necessary to approximate the high-order complex system by low-order simple model. The resulting low-order model is less expensive and more efficient computationally and hence it is used as a replacement to the original high-order model. Over the last few decades, a great deal of attention has been paid to model reduction techniques and many model reduction algorithms
∗ Tel.: +966 2 695 2195; fax: +966 2 695 2686.
E-mail address: [email protected]. 1434-8411/$ - see front matter 䉷 2005 Elsevier GmbH. All rights reserved. doi:10.1016/j.aeue.2005.05.016
have been proposed. Earlier methods were mainly based on retaining the dominant poles in the reduced model as in singular perturbation approximation [1], or on matching some moments of the original system as in Padé approximation techniques [2]. An interesting approach for model reduction is the one in which each component of the state vector has equal information about controllability and observability, i.e., the input-to-state coupling and the state-to-output coupling are weighted equally. This approach is known in the literature as the balanced state–space realization [3–7]. The reducedorder is determined by truncating the system states to those associated with the most significant singular values and the remaining states are rejected. The drawback of balanced
414
Rabah W. Aldhaheri / Int. J. Electron. Commun. (AEÜ) 60 (2006) 413 – 420
realization-based algorithms is that it requires the original system to be completely controllable and observable and consequently, these algorithms are ill-conditioned when the original system is close to being non- minimal. In [8,9], an equivalent to the balanced reduced model is proposed without computing the balancing transformation. Moreover, it was shown in [9] that the balanced reduced model can be found even when the original model is non-minimal. Other interesting approaches for model reduction are the ones based on Hankel-norm optimal approximation [10,11], and impulse response gramians [12]. In digital filter design, the aforementioned model reduction algorithms are used to approximate the high-order filters (often FIR) to reduced-order IIR filters [6,7,11–14]. The major reason for starting with FIR filters is because the requirements of linear phase can always be met in the design. Then the FIR is approximated by IIR filters and the linearity of the phase can be preserved, at least, in the passband region. In these techniques, the maximum magnitude error between the full-order and the reduced order is bounded from above by twice the sum of the singular values associated with the eliminated states. Hence, there is still need for improving the error of the reduced filters. In [13], the singular perturbation technique is shown to improve the error in the low-frequency range, and in [14], the finite and optimal search methods are presented to reduce the error in the passband. In [15], the frequencydomain-balanced structure approach is developed and it was shown that the balanced realization obtained by this approach shares many similar properties with its time-domain counterpart. In this paper, we extend the balanced structure approach and introduce the concept of frequency-domain cross-gramian. Moreover, a modified Lyapunov equation is derived to solve the frequency-domain cross-gramian matrix over any specified frequency range, the band where we would like to reduce the approximation error over it. The eigenspace of this matrix is used to convert the large-order filters to low-order filters. Unlike the method in [15], the proposed method does not solve two Lyapunov equations and avoids computing the balanced transformation. This paper is organized as follows: In Section 2, we review some aspects of balanced discrete systems and give the necessary background to formulate the problem addressed in this paper. In Section 3, the concept of the frequencydomain cross-gramian matrix and how it can be computed is presented. In Section 4, we show how to compute orthogonal bases which span the right eigenspaces associated with the large eigenvalues of the frequency-domain cross-gramian matrix using numerically stable method. The Sylvester equation is used to compute the left eigenspaces. The proposed algorithm is also summarized in this section. Examples of digital filters design are given in Section 5 to illustrate the significance of the proposed method and conclusions are given in Section 6.
2. Preliminaries and problem formulation Consider an nth-order, linear, time-invariant and asymptotically stable discrete-time system described by X(k + 1) = AX(k) + Bu(k), y(k) = CX(k) + Du(k),
(2.1)
where A ∈ Rn×n , B ∈ Rn×1 and C ∈ R1×n . The controllability and observability gramians are defined as follows: WC = WO =
∞ k=0 ∞
Ak BB T (AT )k , (AT )k C T CAk .
(2.2)
k=0
These gramians satisfy the following Lyapunov equations: AW C AT − WC + BB T = 0,
(2.3a)
AT WO A − WO + C T C = 0.
(2.3b)
Fernando and Nicholson [16] demonstrated that the information in the controllability gramian and the observability gramian are essentially contained in the cross-gramian matrix defined by WCO =
∞
Ak BCAk .
(2.4)
k=0
Equivalently, WCO can be computed by solving the Lyapunov equation AW CO A − WCO + BC = 0.
(2.5)
We refer to the controllability, observability and crossgramian defined by (2.2)–(2.5) as the time-domain gramians in contrast to the frequency-domain gramians, which will be defined later in following section. The controllability and the observability gramians are used in many previous works [3–5] to get a balanced system realization {Ab , Bb , Cb } provided that the original system is minimal. This was achieved by finding a transformation matrix, T such that 1 0 . (2.6) T −1 WC T −T = T T WO T = = 0 2 The entries of the diagonal of 1 and 2 are assumed to be in a descending order and 1 = diag(1 , 2 , . . . , r ) and 2 = diag(r+1 , r+2 , . . . , n ).
(2.7)
The balanced realization system is determined by A11 A12 B1 , Bb = T −1 B = Ab = T −1 AT = A21 A22 B2 and Cb = CT = [ C1 C2 ] , (2.8)
Rabah W. Aldhaheri / Int. J. Electron. Commun. (AEÜ) 60 (2006) 413 – 420
where A11 ∈ Rr×r , B1 ∈ Rr×1 and C1 ∈ R1×r . The subsystem {A11 , B1 , C1 , D} is considered as a good reducedorder approximation of the original system if r ?r+1 . This rth-order subsystem is called a direct truncation approximation of the balanced system or shortly balanced truncation (BT) approximation. In [8], an equivalent reduced-order model to the subsystem {A11 , B1 , C1 } is proposed without calculating the balanced transformation matrix, but with calculating a transformation matrix that decomposes the original system into strongly controllable and observable coupled subsystem and weakly coupled subsystem. The strongly coupled is retained and the weakly coupled is discarded. Moreover in [13], the reduced-order singularly perturbed system is calculated to obtain a better error approximation in the low-frequency range. These algorithms are numerically more robust compared with the algorithms based on balanced transformation. Yet, they do not give good approximation behavior over the desired frequency range. This motivates us to look for a new technique to achieve a better approximation in any frequency range that we are interested in. In the next section, we propose a definition for the frequency-domain cross-gramian and we will show how it can be used for getting a reduced-order model which behaves closely enough to the original model in the desired frequency range. The model that we consider in this paper is a model of digital filters of different types.
3. Frequency-domain cross-gramian In this section, the frequency-domain cross-gramian matrix is defined and the modified Lyapunov equation to compute this matrix is derived. This definition is an extension to the frequency-domain controllability and observability gramians given in [15]. Moreover, the reduced-order technique is different and it is not based on balanced realization as we will see later in the following section. Lemma 3.1. For the discrete-time system described by (2.1), the cross-gramian matrix is given by 1 WCO = (ej I − A)−1 BC(e−j I − A)−1 d, (3.1) 2 − where I is an n × n identity matrix. Proof. The proof is a direct application of Parseval’s theorem to Eq. (2.4).
415
is equal to the time-domain cross-gramian represented by Eq. (2.4). Theorem 3.1. The frequency-domain cross-gramian matrix, WF , defined by (3.2) satisfies the Lyapunov equation AW F A − WF = −KBC − BCK ∗ , where
(3.3)
∗
denotes the complex conjugate and 2 1 K= (I − e−j A)−1 d − I, 2 1 4
(3.4)
where = 2 − 1 radians/second. Proof. Adding and subtracting the terms ej WCO A, e−j AW CO and WCO to Eq. (2.5) to obtain − (ej I − A)WCO (e−j I − A) + WCO (e−j I − A)ej + (ej I − A)WCO e−j = BC. (3.5) Multiply Eq. (3.5) from the left by (ej I − A)−1 , and from the right by (e−j I − A)−1 . Then, take the integration 1/2 12 d for both sides to obtain WF = −
WCO + K1 WCO + WCO K1∗ , 2
(3.6)
where
2 1 K1 = (I − e−j A)−1 d. 2 1 Substituting (3.6) into the LHS of (3.3), yields ∗ A − WCO + K1 WCO + WCO K1 A 2 − − WCO + K1 WCO + WCO K1∗ . 2
(3.7)
(3.8)
It is easy to show that AK 1 = K1 A and AK ∗1 = K1∗ A. Therefore, Eq. (3.8) can be rewritten as −
[AW CO A − WCO ] + K1 [AW CO A − WCO ] 2 + [AW CO A − WCO ]K1∗ = − K1 − I BC − BC K1∗ − I 4 4 ∗ = −KBC − BCK = RHS of (3.3).
It can be shown that if 1 = − and 2 = , then the Lyapunov equation of the frequency-domain cross-gramian defined by (3.3) is equivalent to its time-domain counterpart defined by (2.5).
Definition 3.1. The frequency-domain cross-gramian in the frequency range (1 , 2 ) is defined as 2 1 WF = (ej I − A)−1 BC(e−j I − A)−1 d. (3.2) 2 1
4. Model reduction algorithm using the corresponding eigenspaces of the frequency-domain cross-gramian
If the whole frequency range, i.e., (−, ) is selected in the above integration, the frequency-domain cross-gramian
In this section, we consider a transformation that decomposes the cross-gramian matrix, WF (or WCO ) into two
416
Rabah W. Aldhaheri / Int. J. Electron. Commun. (AEÜ) 60 (2006) 413 – 420
blocks, one with large eigenvalues and the other with small eigenvalues. This transformation separates the strongly coupled controllable and observable states from the weakly coupled controllable and observable states based on the large and small eigenvalues of the cross-gramian matrix. This transformation block diagonalizes the matrix WF into two blocks: one with large eigenvalues and the other with small eigenvalues, as follows: S 0 V −1 WF V = S = L , (4.1) 0 SS where SL and SS are r × r and (n − r) × (n − r) matrices, respectively. i = |i (SL )|, i = 1, 2, . . . , r, are the large (in magnitude) eigenvalues of WF , and i = |i (SS )|, i = r + 1, r + 2, . . . , n are the small (in magnitude) eigenvalues of WF . Now, let the matrices V and V −1 be partitioned as UL −1 , (4.2) V = [ VL VS ] and V = US where UL [ VL US
I VS ] = r 0
0 In−r
.
VS ] = [ C 1
C2 ]
and
D = D.
(4.4)
The reduced-order model is defined as {A11 , B 1 , C 1 , D} and the transfer function of the reduced model is defined as Hr (z) = D + C 1 [zI r − A11 ]−1 B 1 .
W˜ 11 P − P W˜ 22 + W˜ 12 = 0.
(4.7)
Notice here, the Sylvester equation is computed for two upper triangular forms, W˜ 11 and W˜ 22 as a byproduct of the ordered real Schur form. This will reduce the computational effort required for computing the matrix P . Ir P Step 4: The matrix block diagonalizes W˜ F in 0 In−r
(4.6) and puts it in the form of (4.1). Thus if we partition Q in (4.6) as Q=[ Q1 Q2 ], where Q1 ∈ Rn×r , then VL =Q1 and UL = QT1 − P QT2 . Step 5: The reduced-order model state–space representation is defined as 1. 2. 3. 4.
A11 = UL AV L B 1 = UL B C 1 = CV L D = D.
(4.3)
The columns of VL and the rows of UL span the right and the left eigenspaces associated with (SL ), i.e., the large eigenvalues of WF . Similarly, the columns of VS and the rows of US span the right and the left eigenspaces associated with (SS ), i.e., the small eigenvalues of WF . Applying the transformation, Eq. (4.2), to the original system (2.1), we obtain UL A11 A12 A[ VL VS ] = A = , US A21 A22 UL B1 B= , B= US B2 C = C[ VL
where the eigenvalues of WF appear in descending order of absolute value along the diagonal of W˜ F , W˜ 11 and W˜ 22 are r × r and (n − r) × (n − r) matrices, respectively. Step 3: Compute the matrix P ∈ Rr×(n−r) by solving the Sylvester equation
(4.5)
In the rest of this section, let us summarize the proposed algorithm to obtain a reduced-order model and in the next section we apply this algorithm to design filter-order reduction in an arbitrary frequency range. Step 1: Given the original state–space representation {A, B, C, D} of the full-order digital filter, determine the frequency range, (1 , 2 ) over which it is desired to approximate the full-order filter by a reduced-order one. Then, compute WF either by solving the Lyapunov equation as in Eq. (3.3), or equivalently, by numerical integration of (3.2). Step 2: Compute the ordered real Schur form of WF [8,17]. W˜ 11 W˜ 12 QT WF Q = W˜ F = , (4.6) 0 W˜ 22
Substitute into (4.5) to obtain the transfer function of the reduced-order filter Hr (z).
5. Numerical examples In this section, we present three examples of filter-order reduction using the proposed method and two of the widely used time-domain-based algorithms: the BT and the optimal Hankel-norm approximation (OHA) methods to illustrate the effectiveness of the proposed method in obtaining lower magnitude error in the desired frequency range. Our concern in the first two examples is the magnitude response and how to get a better error performance. In the third example the phase response is discussed and it is shown that the phase linearity can always be preserved provided that the full-order filter is of linear phase. The design procedure consists of two steps: first we design a filter that meets certain specifications using one of the standard procedures. Second, we approximate the full-order filter with a reducedorder one using the algorithm proposed in Section 4 such that the aforementioned specifications are met, at least, in the chosen frequency band. Example 1. Design a Butterworth lowpass that meets the following specifications: 1, 0 0.45, j |H (e )| = 0, 0.54 . The minimum stopband attenuation is 25 dB. It is found that an 12th-order IIR filter meets the above specifications. The magnitude response, which we will
Rabah W. Aldhaheri / Int. J. Electron. Commun. (AEÜ) 60 (2006) 413 – 420 -3
1.4 4
Full order, n = 12 BT method, r = 7 OHA method, r = 7 Proposed method, r = 7
1.2
x 10
3 2 Magnitude response of error
1 Magnitude response
417
0.8
0.6
0.4
1 0 -1 -2 -3
0.2
BT method, r = 7 OHA method, r = 7 Proposed method, r = 7
-4
0
0
0.1
0.2
0.3
0.4 0.5 0.6 Normalized frequency
0.7
0.8
0.9
-5
1
Fig. 1. Magnitude responses of the original and reduced-order IIR filters.
0.2
0.3 0.4 Normalized frequency
0.5
0.6
0.7
1.4
x 10
BT method, r = 7 OHA method, r = 7 Proposed method, r = 7
4.5
0.1
Fig. 3. Magnitude responses error of the two reduced-order filters.
-3
5
0
Full order, n = 6 BT method, r = 4 OHA method, r = 4 Proposed method, r = 4
1.2
4
1 Magnitude response
Error magnitude
3.5 3 2.5 2
0.8
0.6
0.4
1.5 1
0.2
0.5 0
0 0
0.1
0.2
0.3 0.4 Normalized frequency
0.5
0.6
0.7
0
0.1
0.2
0.3
0.4 0.5 0.6 Normalized frequency
0.7
0.8
0.9
1
Fig. 2. Error magnitude of the reduced-order IIR filters.
Fig. 4. Magnitude responses of the original and reduced-order highpass IIR filters.
consider it here as the full order, is depicted in Fig. 1. By applying the algorithm proposed in Section 4, a 7th-order IIR filter is selected as a good approximate to the full-order filter over the frequency band [−0.55, 0.55]. Moreover, a 7th-order IIR filter is computed using the conventional BT and OHA methods. Fig. 1 shows the magnitude responses of the full order (dashed curves), the proposed reduced order (solid curves), and the BT method (dotted curves), and the OHA (dashed-dotted curves) are plotted over the normalized frequency, where the normalized frequency 1 represents
the Nyquist frequency. As we notice, the magnitude response of the proposed method is much closer to the full order than the reduced order obtained by the BT and OHA methods. Figs. 2 and 3 show error magnitude and magnitude response error over the normalized frequency band [0, 0.55]. These figures show clearly that the error over the selected band for the proposed method is much better than that for the BT and OHA methods.
H (z) =
Example 2. Consider the highpass IIR filter represented by the transfer function
4.7297 × 10−2 (1 − 6z−1 + 15z−2 − 20z−3 + 15z−4 − 6z−5 + z−6 ) . 1 − 0.486z−1 + 1.2052z−2 − 0.042145z−3 + 0.36827z−4 + 0.08425z−5 + 0.04481z−6 The magnitude response of this filter is shown in Fig. 4. Again, we applied the BT and OHA methods, and the
418
Rabah W. Aldhaheri / Int. J. Electron. Commun. (AEÜ) 60 (2006) 413 – 420 0.08
10
BT method, r = 4 OHA method, r = 4 Proposed method, r = 4
FIR filter, n = 45 BT method, r = 15 OHA method, r = 15 Proposed method, r = 15
0
0.04
-10
0.02
-20
Magnitude in dB
Magnitude response error
0.06
0 -0.02
-30 -40 -50
-0.04
-60
-0.06
-70
-0.08 0.3
0.4
0.5
0.6 0.7 Normalized frequency
0.8
0.9
1
Fig. 5. Magnitude response error of the three reduced-order filters.
0
0.1
0.2
0.3
0.4 0.5 0.6 Normalized frequency
0.7
0.8
0.9
1
Fig. 6. Magnitude responses of lowpass FIR and reduced-order IIR filters.
proposed method with selected frequency range [0.35, ] to obtain a 4th-order IIR filters represented, respectively, by HBT (z) =
10−2 (4.7297 − 27.255z−1 + 64.4025z−2 − 71.306z−3 + 38.592z−4 ) , 1 − 0.28218z−1 + 0.80547z−2 + 0.10036z−3 + 0.14751z−4
HOHA (z) =
10−2 (4.7297 − 26.684z−1 + 63.728z−2 − 69.916z−3 + 38.805z−4 ) 1 − 0.26310z−1 + 0.84277z−2 + 0.10496z−3 + 0.19029z−4
and Hr (z) =
10−2 (4.7297 − 27.307z−1 + 72.732z−2 − 81.597z−3 + 56.266z−4 ) . 1 − 0.53423z−1 + 0.83175z−2 − 0.032746z−3 + 0.10221z−4
The magnitude responses of the lower-order filters obtained by the BT, OHA and the proposed methods are shown in Fig. 4. To compare between the three approximations, the magnitude response error of the BT method and OHA, which is based on time-domain model reduction, and the proposed method, which is based on frequency-domain model reduction, are plotted in Fig. 5. It is shown that the proposed method gives better error performance than that of the BT and OHA methods in the passband [0.4,1]. It is worth mentioning here that the order reduction ratio from FIR to IIR filters is much more than the reduction ratio when we approximate the IIR by another IIR filters. This will be addressed in the following example. Example 3. Design a linear-phase lowpass FIR filter with the following specifications: j
|H (e )| =
1, 0,
0 0.3, 0.35 ,
where the maximum passband attenuation is 1 dB and the minimum stopband attenuation is less than 25 dB. Using the Remez algorithm [18], it is found that an FIR of order 45 meets the above requirements. The magnitude response of this filter is shown in Fig. 6. Now, if we apply the proposed algorithm with a selected frequency band [−0.36, 0.36], an IIR filter of order 15 is a good approximation which tracks the FIR responses fairly closely over the selected band. The magnitude responses of the lowerorder filters obtained by the BT, OHA and the proposed method are shown in Fig. 6. A comparison between the three approximations is depicted in Fig. 7. As it is clearly shown in this figure, the proposed method gives a better error performance over the selected passband [0 − 0.36] compared with that of the BT and OHA methods, respectively. Fig. 8 shows the corresponding group delay for both the FIR filter and the three reduced-order approximations. Notice that the phase obtained by the proposed algorithm is linear over the band [0 − 0.36]. This confirms that the proposed algorithm maintains the phase linearity as it is the case with the time-domain-based algorithms.
Rabah W. Aldhaheri / Int. J. Electron. Commun. (AEÜ) 60 (2006) 413 – 420
phase linearity is preserved over the selected band provided that the original filter has a linear phase.
0.025 BT method, r = 15 OHA method, r = 15 Proposed method, r = 15
0.02
419
Magnitude response error
0.015 0.01
References
0.005 0 -0.005 -0.01 -0.015 -0.02 -0.025
0
0.05
0.1
0.15 0.2 0.25 Normalized frequency
0.3
0.35
0.4
Fig. 7. Magnitude responses error of the three approximated IIR filters.
60 FIR filter, n = 45 BT method, r = 15 OHA method, r = 15 Proposed method, r = 15
40
Group delay
20
0
-20
-40
-60
0
0.1
0.2
0.3
0.4 0.5 0.6 Normalized frequency
0.7
0.8
0.9
1
Fig. 8. Group delay responses for the FIR and the three approximated IIR filters.
6. Conclusions In this paper, the frequency-domain cross-gramian, WF is derived and computed by solving a Lyapunov equation (3.3). The left and the right eigenspaces associated with the large eigenvalues of WF are used to transform the full-order to a reduced-order models without computing the balancing transformation which tends to be ill-conditioned, especially when the original model is non-minimal or when it has very nearly uncontrollable or unobservable modes. This technique is applied to several digital filters designs and it is found that the error performance of the filters design by this technique outperforms the time-domain-based algorithms such as the balanced truncation method and the optimal Hankel-norm approximation over the selected frequency range. Moreover, using the proposed technique, the
[1] Kokotovic P, Khalil H, O’Reilly J. Singular perturbation methods in control: analysis and design. New York: Academic Press; 1986. [2] Shamash Y. Model reduction using the Routh stability criterion and Pade approximations. Int J Control 1975;21: 475–84. [3] Moore BC. Principal component analysis in linear systems: controllability, observability and model reduction. IEEE Trans Automat Control 1981;AC-26:17–32. [4] Pernebo L, Silverman LM. Model reduction via balanced state–space representation. IEEE Trans Automat Control 1982;AC-27:382–7. [5] Fernando KV, Nicholson H. On the structure of balanced and other principal representations of SISO systems. IEEE Trans Automat Control 1983;AC-28:228–31. [6] Beliczynski B, Kale I, Cain GD. Approximation of FIR by IIR digital filters: an algorithm based on balanced model reduction. IEEE Trans Signal Process 1992;40:532–42. [7] Kale I, Gryka J, Cain GD, Beliczynski B. FIR filter order reduction: balanced model truncation and Hankel-norm optimal approximation. IEE Proc Vision Image Signal Process 1994;141:168–74. [8] Aldhaheri RW. Model order reduction via real Schur form decomposition. Int J Control 1991;53:709–16. [9] Aldhaheri RW, Bamani A. Balanced realization of nonminimal MIMO discrete time systems. Int J Systems Sci 1994;25:173–82. [10] Glover K. All optimal Hankel norm approximation of linear multivariable systems and their L∞ error bounds. Int J Control 1984;39:1115–93. [11] Chen B-S, Peng S-C, Chiou B-W. IIR filter design via optimal Hankel norm approximation. IEE Proc G Circuits Devices Systems 1992;139:586–90. [12] Sreeram V, Agathoklis P. Design of linear-phase IIR filters via impulse-response gramians. IEEE Trans Signal Process 1992;40:389–94. [13] Aldhaheri RW. Design of linear-phase IIR digital filters using singular perturbational model reduction. IEE Proc Vision Image Signal Process 2000;147:409–14. [14] Li L, Xie L, Yan W-Y, Soh YC. Design of low-order linearphase IIR filters via orthogonal projection. IEEE Trans Signal Process 1999;47:448–57. [15] Wang D, Zilouchian A. Model reduction of discrete linear systems via frequency-domain balanced structure. IEEE Trans Circuits Systems 2000;47:830–7. [16] Fernando KV, Nicholson H. On the cross-gramians for symmetric MIMO systems. IEEE Trans Circuits Systems 1985;CAS32:487–9. [17] Laub AJ. A Schur method for solving algebraic Riccati equations. IEEE Trans Automat Control 1979;AC-24:913–21. [18] Proakis JG, Manolakis DG. Digital signal processing: principles, algorithms, and applications. 3rd ed., Englewood Cliffs, NJ: Prentice-Hall; 1996.
420
Rabah W. Aldhaheri / Int. J. Electron. Commun. (AEÜ) 60 (2006) 413 – 420
Rabah W. Aldhaheri was born in Jeddah, Saudi Arabia in 1953. He received his B.S. degree from Riyadh (currently, King Saud) University in 1976, the M.S. degree from Ohio University, in 1981, and the Ph.D. degree from Michigan State University in 1988, all in Electrical Engineering. From 1976 to 1981, he was with the P.T.T. Ministry, Jeddah, Saudi Arabia as an Electrical Engineer and then as a
Director of the Microwave and Multiplex Stations. In 1981, he joined King Abdulaziz University as a Lecturer in the Department of Electrical and Computer Engineering. He currently is an Associate Professor in the same department. Dr. Aldhaheri has held visiting research scholar positions with Michigan State University during 1994–1995, and Queensland University of Technology (QUT) in Brisbane, Australia in 2000. His research interests include digital signal processing with application to filter design, speaker recognition, face recognition and iris recognition; and digital communications. He is an Associate Editor of the Journal of Applied Sciences.