Hankel-Norm Approximation of IIR by FIR Models - Semantic Scholar

Report 8 Downloads 95 Views
586

IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 55, NO. 2, MARCH 2008

Hankel-Norm Approximation of IIR by FIR Models: A Constructive Method Li Chai, Member, IEEE, Jingxin Zhang, Member, IEEE, Cishen Zhang, Member, IEEE, and Edoardo Mosca, Life Fellow, IEEE

Abstract—This paper presents a constructive method for (sub)optimal finite-impulse response (FIR) approximation of infinite-impulse response (IIR) models. The method minimizes the Hankel norm of the approximation error by using an explicit solution to the norm-preserved dilation problem. It has advantages over the existing methods in that it is a unified method for both single-input single-output and multiple-input multiple-output systems which allows direct tradeoff between the least-squares and Chebyshev error criteria by using a single tuning parameter, and that the approximation algorithm is constructive and only involves algebraic manipulations rather than iteration and convex optimization procedures. The lower and upper bounds on the 2 and Chebyshev norms of the approximation error are derived and are related to the tuning parameter. The problem of approximating noncausal IIR models by causal FIR models is also formulated and solved. The effectiveness and properties of the proposed algorithms are demonstrated by examples. Index Terms—Finite-impulse response (FIR) approximation, Hankel norm, infinite-impulse response (IIR) models, mixed-norm design, norm-preserved dilation.

I. INTRODUCTION

F

INITE-IMPULSE response (FIR) models have advantages in intrinsic stability, simple structure and easy implementation, thus are more preferred to infinite-impulse response (IIR) models in practical applications such as signal processing [1], [2], telecommunications [3] and control systems [4], [5]. However, IIR models often arise from the modeling of systems and signals and the design of controllers and filters [1], [2], [6]. Therefore, effective methods are required to approximate IIR models by FIR models. In general, such an approximation problem can be stated as follows: Given and IIR , which is a stable rational transfer function, find model such that the norm

Manuscript received April 28, 2006; revised April 21, 2007. This work was supported in part by the Australian Research Council’s Discovery funding scheme (project number DP0343057) and in part by National Science Foundation of China under Grant 60304011 and Grant 60672064. This paper was recommended by Associate Editor C. Hadjicostis. L. Chai is with the Engineering Research Center of Metallurgical Automation and Measurement Technology, Ministry of Education, Wuhan University of Science and Technology, Wuhan 430081, China (e-mail:[email protected]). J. Zhang is with the Department of Electrical and Computer Systems Engineering, Monash University, Vic. 3800, Australia (e-mail: jingxin.zhang@eng. monash.edu.au). C. Zhang is with the School of Electrical and Electronic Engineering, Nanyang Technological University, Singapore 639798 (e-mail: ecszhang@ntu. edu.sg). E. Mosca is with the Dipartimento di Sistemi e Informatica, Universityá di Firenze, Florence 50139, Italy. Digital Object Identifier 10.1109/TCSI.2007.913615

of the error transfer function is minimized, where the specific form of the norm depends on the chosen approximation criteria. Similar approximation problem also arises in FIR filter design, where FIR filters are used to approximate the desired frequency responses which are normally not rational transfer functions but some piecewise-constant and brick-wall form specifications of desired frequency responses. Such design problems have been studied extensively in the literature and many novel and effective methods have been obtained [7]–[14]. All these methods are frequency domain methods which use iterative frequency domain interpolation to find the required FIR approximation of the desired frequency response. In principle, these methods can also be used for FIR approximation of arbitrary IIR rational transfer functions with complicated forms of frequency responses. However, when the frequency response is complicated, the interpolation grid must be dense in order to capture the characteristics of the desired frequency response. This can significantly increase computation complexity of iterative frequency domain interpolation and make the frequency domain methods inefficient. Also, all these methods are developed for single-input single-output (SISO) systems, which cannot be used directly to obtain the optimal FIR approximation of the rational transfer functions with multiple-input multiple-output (MIMO) variables. Therefore, these frequency domain methods may not be best for the FIR approximation of arbitrary IIR rational transfer functions with complicated frequency responses. Most of the existing works on FIR approximation of IIR rational transfer functions use time-domain methods. The simplest one is the direct truncation of the impulse response that minimizes the least-squares error criterion, or equivalently the error norm [2]. In [15]–[17], the minimum Chebyshev error error norm was criterion that minimizes the Chebyshev used. In [16], [17], a method called Nehari Shuffle was proposed to minimize the Chebyshev error norm, and the upper and lower bounds on the approximation error were derived. However, the Nehari Shuffle did not provide optimal solutions with respect to Chebyshev norm. A direct Chebyshev norm optimization approach was recently derived in [15] using linear matrix inequalities (LMIs). This approach can achieve nearly optimal Chebyshev error norm. All these methods are developed for MIMO systems, including SISO systems as a special case. But none of these methods is capable of tradeoff approximation. Apart from the above problems, all these time-domain approximation methods require iterative calculations, which are computationally expensive and may suffer from slow convergence (see Section V for demonstrative examples). To our best knowledge, there has been no time-domain direct tradeoff

1549-8328/$25.00 © 2008 IEEE

CHAI et al.: HANKEL-NORM APPROXIMATION OF IIR BY FIR MODELS: A CONSTRUCTIVE METHOD

method for FIR approximation of IIR models, and there has been no time-domain approximation algorithm that does not involve iterative calculations. To overcome the difficulties of time-domain approximation methods discussed above, this paper considers the Hankel norm minimization of the approximation error. Hankel-norm approximation has been extensively studied in model reduction after the remarkable work of Glover [18], [19]. However, the problem considered in this paper is different from that of [18] where a lower order IIR model is constructed for a given high-order IIR model and rather complicated balanced realization and truncation are required. The method proposed in this paper has the following advantages. • It is developed for MIMO systems, including SISO systems as a special case. • It allows direct tradeoff between the least-squares criterion and Chebyshev criterion by using a single tuning parameter. • The design algorithm is constructive and only involves algebraic manipulations rather than iteration and convex optimization procedures (such as LMIs). • There is no need to carry out balanced realization and truncation as in [17]. • Lower and upper bounds on the and Chebyshev norms of the approximation error are provided and are related to the tuning parameter. The remainder of the paper is as follows. Section II provides some necessary background on Hankel operators and norm-preserved dilations. The basic approximation algorithm is developed in Section III and its extension to noncausal systems is given in Section IV. Section V discusses how to incorporate frequency weighting in the approximation algorithm. Computation examples are given in Section VI to show the effectiveness and properties of the proposed method. Finally the paper is concluded in Section VII. II. PRELIMINARY This section introduces notations and preliminaries used in the sequel. Let and denote the real and complex denote its numbers, respectively. For a matrix , let its eigenvalue, and complex conjugate transpose, its th singular value. Denote the Frobenius norm of as , and the spectrum norm of as , where denotes the largest eigenvalue to denote of . For a positive definite matrix , we use and its Hermitian square root which satisfies .

587

and can be computed from the It is well known that following Lyapunov equations, respectively: (1) (2) The realization is minimal if and are nonsingular [20]. with For causal and stable minimal state space realization , where is the impulse response of at time , , , and , the norm of , denoted by , is given by (3) The

norm (or Chebyshev norm) of , is given by

, denoted by

(4) It is assumed that the right hand sides of (3) and (4) are all well-defined for simplicity. For more rigorous definitions, , please refer to [20]. The Hankel operator of , denoted by is defined as

.. .

.. .

.. .

..

.

The Hankel singular values of , denoted by , , are the th singular values of . The Hankel norm is defined to be the greatest singular of , denoted by value of , i.e., . It is worth noting that and hence are unrelated to the constant term of . The following can be used to compute the Hankel norm of a transfer function, see [18], [20], [21] for details. the following holds. Lemma 1: For the above

where and are the controllability and observability Gramians, respectively. B. Norm-Preserved Dilations Consider the block matrix

, where

,

,

and

A. Spaces, Norms, and Hankel Operators

are matrices of compatible dimensions, and denote

, Definition 1: Given a causal transfer function is called a state space realization of if , where , , and . Definition 2: For a stable transfer function with the state , the controllaspace realization bility and observability Gramians, denoted by and , are defined by and .

The norm-preserved dilation problem is to find such that is minimized for the given matrices , , and . Denote (5)

588

IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 55, NO. 2, MARCH 2008

The following results [22] play a very important role in our development. Lemma 2: The in (5) is given by

Moreover, for a given , the solution set can be characterized by

such that (6)

is an arbitrary contraction with where are contractions satisfying

, and

, define

Theorem 1: For . Then we have

where and are solutions of the Lyapunov (1) and (2), respectively. can be written in the state space Proof: Note that , where form as

and (7) (8)

The following lemma gives a more explicit formula when . and . Then one Lemma 3: Assume that is given by solution of such that

Thus, for

, the Lyapunov (1) and (2) become (12) (13)

Writing

as

and substituting it into (12)

yield

(9) The norm-preserved dilation problem is solved independently by Parrot and Davis et al. For more details, please refer to [22].

From (1)

III. HANKEL-NORM FIR APPROXIMATION A. Development of the Algorithm This section develops the algorithm for Hankel norm FIR approximation. The problem considered is as follows. Given an , find an IIR model th-order FIR model

It then follows that

Similarly, writing

as

and substituting it

into (13) yield with , that minimizes , the Hankel norm of , where is the approximation error. The reason to consider the delayed error is that the Hankel norm of a transfer function is unrelated to its constant term, and that the delay operator preserves the and Chebyshev norms. The relation of the Hankel norm, norm and Chebyshev norm for and are given in the following lemma. Lemma 4: Let and with . Then (10)

It then follows from (2) that

Using direct matrix manipulation, it is easy to show that

(11) where

is the th Hankel singular value of and . Proof: The equalities in (10) follow directly from the definition of and Chebyshev norms. The first two inequalities in (11) are shown by [23, Th. 4.2], and the last inequality is shown in [17], [20]. Lemma 4 shows that the Hankel norm can be seen as a tradeoff between the norm and Chebyshev norm. The following theorem is important to the development of new approximation algorithm.

Then it follows from Lemma 1 that

CHAI et al.: HANKEL-NORM APPROXIMATION OF IIR BY FIR MODELS: A CONSTRUCTIVE METHOD

Theorem 2: Given a transfer function with , define . Then the following holds. for any . i) ii) There exist ’s such that can be characterized by

for a matrix

, and all such

’s

589

By Lemma 3, the satisfying the above equation gives rise to . This proves part iii). Remark 1: Theorem 2 is instrumental to the development of Hankel norm approximation algorithm. To develop the algorithm using the theorem, we recall the following well known fact can [17]: a causal transfer function be written in the form , where

(14) where

iii) For any one of such

and

and

(20)

are contractions satisfying

, there exist ’s such that ’s is given by

, and

(15) Proof: It follows from Theorem 1 that

, and with being the first impulse responses of is a strictly proper transfer function. In the direct truncation as given in (20) is used as the th-order method, FIR approximation of , and is the truncation error that has the minimum norm among all th-order FIR approximations. The theorem , an FIR filter below will show that from can be constructed such that for a prescribed satisfying , the delayed approximation error defined below attains .

(21) Note that

(16)

(17)

be Theorem 3: Let a stable and causal transfer function with , and let with . Then the following holds. for any . i) ii) For any , a particular can be constructed explicitly such that . Proof: Denote and

By Lemma 1, (16) and (17) (22) for gives (18) It then follows from Lemma 2 that

. Using part i) of Theorem 2 recursively

(23) for any . This proves part i). To prove part ii), a constructive approach will be used to show can be computed if is obtained. Now assume that that a state-space realization for is given by

Moreover, there exists such that . Substituting , and into the formula (6), it is easy to obtain (14) after some direct algebraic manipulations. This proves (i) and (ii). Obviously . Substituting , and into (9) gives

(19)

(24) Then the controllability and observability Gramians and can be computed from (1) and (2), respectively. Since , it then follows from part iii) of Theorem 2 that there exists such that

590

where given by

IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 55, NO. 2, MARCH 2008

is defined by (22). Moreover, one

is

(25) Following the same line as the proof of Theorem 1, it is easy to is given by check that a state space realization for (26)

5) Obtain

by the following equation

(31) 6) Obtain a state space realization of from (26) and . and obtain 7) Set to and repeat the steps 5) and 6) to find , for . 8) The optimal Hankel-norm approximation is then given by

where (32)

and

The controllability and observability Gramians and for the state space realization (26) are as shown in (27) and (28) at the bottom of the page. The proof is then completed by Theorem 2 by noting that we can now compute again. Remark 2: Substituting (26)–(28) to (25), the following formula with lower computational load can be obtained after some straightforward operations.

Remark 3: Algorithm 1 only provides a suboptimal solu. Setting and using (14), the option with can also be obtained, although it timal solution with is not explicit. Actually, the optimal solution can be obtained from Algorithm 1 by simply using the pseudo-inverse of if it is singular. B. Properties of the Algorithm The following Corollary gives lower and upper bounds on the and Chebyshev norms of the approximation error of Algorithm 1, and reveals the relationship between these bounds and the parameter in the algorithm. Corollary 1: For with , let be the FIR approximation of obtained by Algorithm 1. Then the following holds for the . approximation error

(29)

(33)

and . The formula holds for The proof of Theorem 3 shows how the direct truncation can be modified to obtain the th-order (sub)optimal Hankel norm FIR approximation of a given IIR filter . The procedure is summarized in the algorithm below. Algorithm 1

(34)

(35) where

1) Decompose

in the form with

as given in

(20). 2) Set

, where and . 3) Obtain and by solving the Lyapunov (1) and (2). 4) Compute using any of the following equations, and specify a . ,

is the th Hankel singular value of , and is a finite constant. Proof: The inequalities in (33) and (34) follow directly from Lemma 4 and Theorem 3. To prove (35), consider (32) and (21). It is easy to see that and . Denote . The first two inequalities in (35) then follow from Lemma 4, Theorem 3 and the fact that . From Lemma 1 and (23), for ,

(30)

(27) (28)

CHAI et al.: HANKEL-NORM APPROXIMATION OF IIR BY FIR MODELS: A CONSTRUCTIVE METHOD

Denote

. Then we have . It then follows from

(29) that

(36) where

is a constant satisfying .

By

Lemma 4, . Thus, all the terms in

for

are . Because (36) holds for all , the last inequality in (35) then follows . from taking Corollary 1 shows that the upper bound on the Chebyshev norm of approximation error is . Actually a tighter upper bound can be attained simply by another choice of . The result is as follows, see [18], [20] for details. Corollary 2: For , let be chosen as in Algorithm 1. If is chosen such that is minimized, then we have bounded, and

(37) Remark 4: Corollaries 1 and 2 have revealed a very important property of the parameter in Algorithm 1. In fact, is a direct tuning parameter for the tradeoff between the least squares and Chebyshev criteria. This can be easily observed from (33)–(35) , the smaller the and (37). Generally speaking, for all is, the smaller the is, and vice versa. On the other hand, is, and vice versa. Furthe larger the is, the smaller the thermore, if , then , the minimum error norm attained by direct truncation. will be further demonstrated This tradeoff function of through computation examples in Section VI. IV. APPROXIMATION OF NONCAUSAL SYSTEMS There are situations where causal filter (transfer function) approximation of noncausal filters are required [17], [24]. This section considers the problem of approximating a noncausal by a causal FIR filter. First we consider the anticausal case. For an anticausal IIR filter , denote . It can be the reverse operator as readily shown that the reverse operator preserves norms. Hence, the Hankel, and Chebyshev norms and the th Hankel singular value of are given, respectively, by , , , and . The problem is stated as follows: given an anticausal IIR filter , find an FIR filter

591

that minimizes , where is the Hankel norm of the forward shifted approximation error . The reason to consider the forward shifted error is the same as that for causal case. The theorem below shows that this problem can be easily solved by converting the anticausal approximation problem to an equivalent causal approximation problem. Theorem 4: Given an anticausal IIR filter with , there exits a causal FIR approxsuch that , imation filter is the Hankel norm of the approximation error where . Moreover, such an can be obtained by constructing an FIR filter that achieves and by taking , is the Hankel norm of the reversed approximawhere tion error . Proof: By definitions, , and . Let , for . Then can be written in the form (21), namely, . Note is causal and strictly that . It then folproper, and that lows from Theorem 3 that a and hence an can be constructed to achieve . Because , the can be reversed to obtain constructed . The thus . obtained is causal and FIR, and yields This proves the theorem. Theorem 4 shows that Algorithm 1 can be used to find the (sub)optimal Hankel-norm FIR approximant for a given anticausal IIR filter. This is summarized in Algorithm 2 below. Algorithm 2 1) Obtain a state space realization of . , where 2) Set and . . 3) Use Algorithm 1 to find 4) Obtain the approximation as . Corollary 3: For by Algorithm 2, let Then

, let and

,

be obtained .

Proof: The results are obvious from Corollary 1 by noting that the reverse operator preserves the norms. So far, we have developed algorithms to obtain the Hankelnorm FIR approximation of causal and anticausal IIR filters, re-

592

IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 55, NO. 2, MARCH 2008

spectively. For a general noncausal IIR filter, we have the following algorithm using the treatment in [17], [24].

and

Define

Algorithm 3

(39)

into its causal 1) Decompose the noncausal and anticausal (both stable) components: . 2) Use algorithm 1 to get the Hankel-norm FIR to . approximation 3) Use algorithm 2 to get the Hankel-norm FIR to . approximation 4) Obtain the approximation as Corollary 4: For tained by Algorithm 3, let , Then

, let

with

can be further written as (40)

is a strictly proper transfer function with and that (40) is in the same form as (21). By ,a can be conTheorem 3, for any . Therefore, Algorithm 1 structed such that is invertible, the can be used to obtain such a . Because s that give can be solved by Notice that

. be ob-

, and

. Then

.

(41)

Proof: Note that

The results then follow directly from Corollary 1 and Corollary 3.

V. FREQUENCY WEIGHTED HANKEL-NORM FIR APPROXIMATION The approximation method presented in Section III can be extended to incorporate frequency weighting although not straightforward. This section discusses how this can be done and the associated difficulties and possible solutions. To simplify the presentation, only the input weighting will be considered here. The results carryover easily to the output weighting and the combination of input and output weightings. be an invertible frequency weighting function of Let the form with . Define the frequency weighted approximation error and the delayed frequency weighted approximation . Using (21), error can be written as

The Hankel norm FIR approximation now is to find minimizes .

(38) s that

is the transfer function defined in (39). Substituting where the solved s into (32) gives the required FIR approximation . needs to As seen from above, the weighting function be FIR, and its inverse is required in (41) to find the optimal FIR approximation. These impose some constraints on the de. The order of will impact both the approxisign of mation error and the computation complexity. A careful design is therefore required to tradeoff between these two factors. The extraction of s from (41) relies on the cancellation of the (transmission) zeros of with the poles of . If these zeros and poles are not exactly same due to computation errors, the precise cancelation may not be possible. For SISO systems , an approximated solution of s for such a case can th-order quotient of the long be obtained by taking the division. For MIMO systems, the approximated solution is yet to be found. All these issues are currently being investigated. VI. EXAMPLES This section presents examples of FIR approximation of IIR models to illustrate the proposed algorithms. Example 1: Consider the following sixth-order IIR filter:

This is the model of spindle vibration we obtained at a hot steel rolling mill for prediction and reduction of mechanical failure [25]. The model is nonminimum phase and has a pole very close to the unit circle in the -plane. Hence, it is prone to numerical error and not suitable for digital signal processing (DSP) implementation. To overcome this implementation difficulty, an FIR approximation is required.

CHAI et al.: HANKEL-NORM APPROXIMATION OF IIR BY FIR MODELS: A CONSTRUCTIVE METHOD

Fig. 1. Example 1: Comparison of frequency response for m = 12.

As shown in Fig. 1 (original), the frequency response of spikes at about , 1.45, 2.02. These spikes, particu, 1.45 cause mechanical damage to larly those two at the spindle [25]. Thus, for this particular application, we need an FIR approximation that captures these two spikes. Now we with the use Algorithm 1 to find FIR approximations of and , respectively. lengths Figs. 1 and 2 compare the frequency response of the origwith those of the FIR approximations obtained by inal Algorithm 1, direct truncation of impulse response, and direct norm optimization [15]. As seen from the figChebyshev ures, when , the FIR approximation tends to smooth optimization attempts to better out all the spikes; while the , 1.45, it generates large deviations capture the spikes at . Tuning to over (note that differs for different , see Remark 1 for details), , 1.45 than Algorithm 1 better captures the spikes at direct truncation and in the same time generates much less dethan norm optimization. viations over When , all the three methods achieve much better ap, Algorithm 1 (with proximation. For the first spike at ) and norm optimization attain almost the same approximation errors, while direct truncation attains a larger error. For the second and third spikes at , 2.02, Algorithm 1 and direct truncation attain almost norm optimization the same approximation errors, whereas attains a much larger error. Apparently, Algorithm 1 achieves the best of both worlds. Fig. 3 compares the and Chebyshev norms of approximation errors at different length attained by Algorithm 1, direct norm optimization. Again, the parameter truncation and in Algorithm 1 is chosen as . As can be seen from the figure, for all , the error norms of these three methods are all below their respective Chebyshev error norms, and the (Chebyshev) error norm of Algorithm 1 is above (below) that of direct truncation and below (above) that of norm optimization. These agree with the analysis of Corollary 1, and demonstrate that Algorithm 1 truly provides a tradeoff between approximation criteria. The figure also the Chebyshev and

593

Fig. 2. Example 1: Comparison of frequency response for m = 24.

Fig. 3. Example 1: l and Chebyshev error norms attained with different m.

demonstrates that the difference of these three methods are significant when the length is smaller. As increases, the difference becomes less significant. When is sufficiently large, these three methods produce the same FIR approximation. Fig. 4 shows the and Chebyshev error norms of the approximation attained by Algorithm 1 with different values, and compares these norms with those attained by direct truncanorm optimization. As can be seen from the figure, tion and , the (Chebyshev) error norm for all satisfying norm opof Algorithm 1 is always below (above) that of timization and above (below) that of direct truncation. When increases from its minimum feasible value , the Chebyshev error norm of Algorithm 1 increases from its minerror norm decreases from its maximum value, while its imum value and approaches the minimum error norm attained by direct truncation as . These verifies the analysis of Corollary 1 and demonstrate that the parameter in Algorithm 1 is truly a direct tuning parameter for the tradeoff between the Chebyshev and least squares approximation criteria.

594

IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 55, NO. 2, MARCH 2008

Fig. 5. Example 2: Comparison of frequency response of three different methods. Fig. 4. Example 1: l and Chebyshev error norms attained with different .

Example 2: This example considers FIR approximation of an 8th-order Chebyshev filter generated by the MATLAB command cheby1(8,0.5,0.3). It is presented here to demonstrate the behavior of Algorithm 1 in approximation of an IIR model with close to brick-wall frequency response, but not to show or suggest the indirect FIR filter design using Algorithm 1. This example is also used in [15], [17]. However, the IIR filters given there are unstable due to the roundoff of filter coefficients. This also shows why IIR filter is not preferred in practical use. Fig. 5 compares the frequency responses of the original Chebyshev filter and the 32-tap FIR filters obtained by Algo), direct truncation and norm rithm 1 (with optimization. As can be seen from the figure, direct truncation gives the best approximation within the pass and stop bands, but its transit band approximation is rather poor, particularly on the pass band boundary. Contrary to direct truncation, norm optimization gives the best approximation on pass band boundary, but its pass and stop band approximation errors are about 20%, which makes it almost useless. Compared with these two methods, Algorithm 1 gives a better approximation than direct truncation on the pass band boundary and a much norm optimization within the better approximation than pass and stop bands. So it has the advantages of both methods. Fig. 6 shows the effects of on the frequency responses of the 32-tap FIR filter obtained by Algorithm 1. As can be seen from the figure, the larger the , the closer the frequency response to that of the FIR filter obtained by direct truncation. Comparing the figure with Fig. 5, it can also be seen that the smaller the , the closer the frequency response to that of the FIR filter obnorm optimization. These again verify the analtained by ysis of Corollary 1 and demonstrate the tradeoff function of the parameter between the Chebyshev and least squares criteria. Fig. 7 compares the and Chebyshev norms of approximation errors at different length attained by Algorithm 1, direct norm optimization. The parameter in Altruncation and gorithm 1 is chosen as . Similar to that of Example 1 shown in Fig. 3, the error norms of these three methods are all below their respective Chebyshev error norms, and the (Chebyshev) error norm of Algorithm 1 is above (below) that of direct truncation and below (above) that of norm op-

Fig. 6. Example 2: Comparison of frequency response attained with different

.

timization. These again agree with the analysis of Corollary 1 and demonstrate that Algorithm 1 truly provides a tradeoff between the Chebyshev and approximation criteria. The figure also reconfirms that the and Chebyshev error norms of these three methods all decrease as increases and all vanish when , and that the difference of these three methods is significant when the length is small and becomes less significant as increases. Fig. 8 shows the and Chebyshev error norms of the approximation attained by Algorithm 1 with different values, and compares these norms with those attained by direct truncanorm optimization. Comparing Fig. 8 with Fig. 4, tion and it can be seen that Algorithm 1 exhibits exactly the same property as in Example 1, particularly, the effects of on the and Chebyshev error norms are exactly the same. This further verifies the analysis of Corollary 1 and demonstrates that Algorithm 1 provides a tradeoff between the Chebyshev and approximation criteria, and that is the direct tuning parameter for the tradeoff between the Chebyshev and least squares approximation criteria. Example 3: This example demonstrates the behavior of Algorithm 1 in approximation of a high-order IIR model generated by MATLAB command butter(30,0.3), and com-

CHAI et al.: HANKEL-NORM APPROXIMATION OF IIR BY FIR MODELS: A CONSTRUCTIVE METHOD

Fig. 7. Example 2: l and Chebyshev error norms attained with different m.

595

Fig. 9. Example 3: Frequency responses of the 30th-order Butterworth and 128-tap FIR approximation filters.

Fig. 8. Example 2: l and Chebyshev error norms attained with different .

pares the computation time of Algorithm 1 with that of LMI norm optimization. A 128-tap FIR filter was combased puted with Algorithm 1 to approximate this 33th-order IIR filter. The frequency responses of the original and computed FIR filters are given in Fig. 9, and the comparison of computation times are given in Remark 6. Example 4: Consider the following IIR transfer matrix:

This is the model of an industrial dry process rotary cement kiln given in [26]. The original model data is given in state space form and is converted to the above transfer matrix. An with and is FIR model computed using Algorithm 1. The singular value frequency reand are compared in Fig. 10. The reason sponses of for such a comparison is that for MIMO systems the singular value frequency responses are the principal gains between inputs and outputs, which are the MIMO system extension of

Fig. 10. Example 4: Singular value frequency responses of G(z ) and F (z ).

SISO system frequency response and are widely used in the analysis and design of MIMO systems [20]. Example 5: This example demonstrates FIR approximation of a two channel IIR synthesis filter bank. The IIR synthesis filter bank is the Example of [27] designed by mixed norm optimization of its polyphase transfer matrix, hence its FIR approximation needs to approximate the principal gains (singular values) of the polyphase transfer matrix of the filter bank. of the synthesis The polyphase transfer matrix filter bank has Smith–McMillan degree of 20. The corresponding synthesis filters are a pair of IIR filters with the order of 41. Using Algois approximated by an FIR transfer matrix rithm 1, the with for and , respectively. The corresponding FIR synthesis filters are 15 taps for and 17 taps for . The singular value frequency responses and the the frequency responses of the FIR synthesis of

596

IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 55, NO. 2, MARCH 2008

TABLE I COMPARISON RESULTS OF EXAMPLES 1–5

Fig. 11. Example 5: Singular value frequency responses of R (z) with m = 9 (dotted).

R(z) (solid) and

norm while pushing the error norm to a rather large value. The gap between the error norms and the gap between the Chebyshev error norms of these two methods are rather large. When is close to its minimum feasible value , Algorithm 1 is similar norm optimization. Its Chebyshev error norm is close to to the minimum value attainable by norm optimization, and its error norm is below that of norm optimization but quite far away from the minimum value attainable by direct truncation. However, with an increased , Algorithm 1 can attain a significantly reduced error norm without much increase of Chebyshev error norm. This is a very desirable property. Because Algorithm 1 involves only algebraic manipulations, it is much faster than the existing time-domain methods [15]–[17] which involve iteration and convex optimization. Tests have been performed using MATLAB on a Toshiba Tecra Notebook computer with a 1.7GHz Pentium(R)M processor to compare the computation time of Algorithm 1 with that of LMI norm optimization [15]. For Example 2, Algorithm based 1 used 0.3710 second to compute the 32-tap FIR approximation norm optimization took a few filter, while LMI based minutes. For Example 3, it took 5.4980 seconds for Algorithm 1 to compute the 128-tap FIR filter, whereas LMI based norm optimization could not obtain any result after hours and had to be abolished. Table I summaries the comparison results of Examples 1–5. VII. CONCLUSION AND DISCUSSION

Fig. 12. Example 5: Frequency responses of the synthesis filters.

filters and are plotted in Figs. 11 and 12, respectively. As can be seen from Figs. 11 and 12, the 17-tap FIR approximations of Algorithm 1 rather precisely capture the frequency characteristics of the original IIR transfer matrices. Compared with the original IIR transfer matrices, the FIR approximation matrices have much less complexity and are numerically more robust due to their intrinsic stability. Remark 5: As can be seen from Fig. 8 and Fig. 4, direct truncation minimizes the error norms while pushing the Chebyshev error norm to a rather large value; opposite to direct trunnorm optimization minimizes the Chebyshev error cation,

A time-domain constructive method is presented to obtain the optimal FIR Hankel norm approximation for a given IIR rational transfer function. The method has advantages over the existing methods in that it is a unified method for both SISO and MIMO systems which allows direct tradeoff between the least-squares and Chebyshev error criteria by using a single tuning parameter, and that the design algorithm is constructive and only involves algebraic manipulations rather than iteration and convex optimization procedures. Hence, the proposed method is practically more useful, efficient and reliable than the existing time-domain methods. The lower and upper bounds on the and Chebyshev norms of the approximation error are derived and their relations to the tuning parameter are established. The problem of approximating noncausal IIR filters by causal FIR filters is also formulated and solved. The effectiveness and properties of the proposed algorithms are demonstrated by examples. The tuning parameter in Algorithm 1 controls directly the tradeoff between the least-squares and Chebyshev error criteria. But it is not a direct measure of either error. Hence, the control is not rigorous. As shown in Section V, the presented method can be extended to incorporate frequency weighting. However, the full solution is yet to be developed. These are weaknesses of the presented method and subjects of our current research.

CHAI et al.: HANKEL-NORM APPROXIMATION OF IIR BY FIR MODELS: A CONSTRUCTIVE METHOD

REFERENCES [1] S. K. Mitra, Digital Signal Processing: A Computer-Based Approach, 2nd ed. New York: McGraw-Hill, 2001. [2] A. V. Oppenheim and R. W. Schafer, Discrete-Time Signal Processing, 2nd ed. Englewood Cliffs, NJ: Prentice-Hall, 1989. [3] M. J. Juntti and B. Aazhang, “Finite memory-length linear multiuser detection for asynchronous CDMA communication,” IEEE Trans. Commun., vol. 45, pp. 611–622, May 1997. [4] C. E. Garcia, D. M. Prett, and M. Morari, “Model predictive control—A survey,” Automatica, vol. 25, pp. 335–348, 1989. [5] W. H. Kwon, “Why not use FIR filters for control?: FIR filters for state space models,” in Proc. 4th Asian Contr. Conf., Singapore, 2002, pp. 6–17. [6] L. Ljung, System Identification, Theory for the Users, 2nd ed. Englewood Cliffs, NJ: Prentice-Hall, 1999. [7] J. W. Adams, “FIR digital filters with least-squares stopbands subject to peak-gain constraints,” IEEE Trans. Circuits Syst., vol. 39, no. 4, pp. 376–388, Apr. 1991. [8] C. S. Burrus and J. Barreto, “Least -power error design of FIR filters,” in Proc. IEEE Int. Symp. Circuits Syst., San Diego, CA, 1992, pp. 545–548. [9] L. J. Karam and J. H. McClellan, “Chebyshev digital FIR filter design,” Signal Process., vol. 76, pp. 17–36, 1999. [10] W. Lertniphonphun and J. McClellan, “Unified design algorithm for complex FIR and IIR filters,” in Proc. ICASSP, May 2001, vol. 6, pp. 3801–3804. [11] C.-Y. Tseng, “An efficient implementation of Lawson’s algorithm with application to complex Chebyshev FIR filter design,” IEEE Trans. Circuits Syst. II, Analog Digit. Signal Process., vol. 42, no. 4, pp. 245–260, Apr. 1991. [12] X. Chen and T. W. Parks, “Design of FIR filters in the complex domain,” IEEE Trans. Acoust., Speech, Signal Process., vol. 35, no. 2, pp. 144–153, Feb. 1987. [13] M. Z. Komodromos, S. F. Russell, and P. T. P. Tang, “Design of FIR filters with complex desired frequency response using a generalized Remez algorithm,” IEEE Trans. Circuits Syst. II, Analog Digit. Signal Process., vol. 42, no. 4, pp. 274–278, Apr. 1995. [14] Y. C. Lim, J. H. Lee, C. K. Chen, and R. H. Yang, “A weighted least squares algorithm for quasi-equiripple FIR and IIR digital filter design,” IEEE Trans. Signal Process., vol. 40, no. 3, pp. 551–558, Mar. 1992. [15] Y. Yamamoto, B. D.O. Anderson, M. Nagahara, and Y. Koyanagi, “Optimizing FIR approximation for discrete-time IIR filters,” IEEE Signal Process. Lett., vol. 10, no. 9, pp. 273–276, Sep. 2003. [16] P. J. Kootsookos, R. R. Bitmead, and M. Green, “The Nehari shuffle: FIR(q) filter design with guaranteed error bounds,” IEEE Trans. Signal Process., vol. 40, no. 8, pp. 1876–1883, Aug. 1992. [17] P. J. Kootsookos and R. R. Bitmead, “The Nehari shuffle and Minimax FIR filter design,” in Control and Dynamic Systems. Orlando, FL: Academic, 1994, vol. 64, pp. 239–298. [18] K. Glover, “All optimal Hankel-norm approximations of linear multivariable systems and their -error bounds,” Int. J. Contr., vol. 39, pp. 1115–1193, Jun. 1984. [19] G. Obinata and B. D. O. Aderson, Model Reduction for Control System Design. Berlin, Germany: Spring-Verlag, 2001. [20] K. Zhou, J. C. Doyle, and K. Glover, Robust Optimal Control. Englwood Cliffs, NJ: Prentice-Hall, 1996. [21] C. Foias, A. E. Frazho, I. Gohberg, and M. A. Kaashoek, Metric Constrained Interpolation, Commutant Lifting and Systems. Berlin, Germany, Birkhäuser, 1998. [22] C. Davis, W. M. Kahan, and H. F. Weinberger, “Norm-preserving dilitions and their applications to optimal error bounds,” SIAM J. Numer. Anal., vol. 19, no. 3, pp. 445–469, Jun. 1982. [23] C. K. Chui and G. Chen, Discrete Optimization: With Application in Signal Processing and Control Systems, 2nd ed. New York: Springer-Verlag, 1997. [24] G. Gu and B. A. Shenoi, “A novel approach to the synthesis of recursive digital filters with linear phase,” IEEE Trans. Circuits Syst., vol. 38, no. 6, pp. 602–612, Jun. 1991.

p

597

[25] K. Zhang, J. Zhang, and S. Nahavandi, “Mechanical failure reduction of hot steel rolling mills via modeling and control,” in TProc. 4th Asian Contr. Conf., Singapore, 2002, pp. 1630–1635. [26] T. Westerlund, “A digital quality control system for an industrial dry process rotary cement kiln,” IEEE Trans. Autom. Contr., vol. AC-26, no. 4, pp. 885–890, Aug. 1981. [27] J. Zhang, R. Yang, C. Zhang, and E. Mosca, “Design of IIR multirate filter banks subject to subband noises,” in Proc 5th Asian Contr. Conf., Melbourne, Australia, Jul. 2004, pp. 1333–1338.

Li Chai (S’00–M’03) received the B.S. degree in applied mathematics and the M. S. degree in control science and engineering, both from Zhejiang University, China, in 1994 and 1997, respectively, and the Ph.D. degree in electrical engineering from Hong Kong University of Science and Technology, Hong Kong, in 2002. In September 2002, he joined Hangzhou Dianzi University, Hangzhou, China. He worked as a Postdoctoral Research Fellow at the Monash University, Clayton, Australia, from May 2004 to June 2006. In December 2007, he joined Wuhan University of Science and Technology, Wuhan, China, where he is currently a Xiang Tao Professor. His research interests include multirate signal processing, wavelets, system identification, and robust control.

Jingxin Zhang (M’02) received the M.Eng. and Ph.D. degrees in electrical engineering, from Northeastern University, Shenyang, China, in 1983 and 1988, respectively. Between 1988 and 1992, he was with Northeastern University as Associate Professor. Since 1989, he has held research positions in University of Florence, Italy, University of Melbourne, and Cooperative Research Center for Sensor Signal and Information Processing, Australia, and Senior Lecturer position in University of South Australia, and Deakin University, Australia. He is currently with the Department of Electrical and Computer Systems Engineering, Monash University, Clayton, Australia. He is the author and coauthor of many research papers in diverse research areas such as adaptive and predictive control, time varying systems, robust filtering, multirate signal processing and medical imaging. His current research interests are in control and signal processing and their applications to industrial and medical systems. Dr. Zhang is the recipient of 1989 Fok Ying Tong Educational Foundation (Hong Kong) for the outstanding Young Faculty Members in China and 1992 China National Education Committee Award for the Advancement of Science and Technology.

l

H

Cishen Zhang (S’87–M’89) received the B.Eng. degree from Tsinghua University, Beijing, China, and the Ph.D. degree in electrical engineering from Newcastle University, Newcastle, NSW, Australia, in 1982 and 1990, respectively. Between 1971 and 1978, he was an Electrician with Changxindian (February Seven) Locomotive Manufactory, Beijing, China. He carried out research work on control systems at Delft University of Technology, Delft, The Netherlands, from 1983 to 1985. After his Ph.D. study from 1986 to 1989 at Newcastle University, he was with the Department of Electrical and Electronic Engineering at the University of Melbourne, Australia as a Lecturer, Senior Lecturer, Associate Professor, and Reader till October 2002. He is currently with the School of Electrical and Electronic Engineering and School of Chemical and Biomedical Engineering at Nanyang Technological University, Singapore. His research interests include signal processing, medical imaging and control.

598

IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS—I: REGULAR PAPERS, VOL. 55, NO. 2, MARCH 2008

Edoardo Mosca (S’63–M’65–SM’95–F’97–LF’06) received the Dr. Eng. degree in electronics engineering from the University of Rome “La Sapienza,” Rome, Italy, in 1963. He then spent four years in the aerospace industry where he worked on the research and development of advanced radar systems: particularly, optimal signal synthesis and processing, and phased-array radar systems. Thereafter, from 1968 to 1972, he held academic positions at the University of Michigan, Ann Arbor, MI, and McMaster University, Hamilton, ON, Canada. Since 1972, he has been with the Engineering Faculty, University of Florence, Florence, Italy: from 1972 to 1974 as an Associate Professor, and since 1975, as a full Professor of Control Engineering. In the latter capacity, he founded in 1981 the Department of Systems and Computer Science and Engineering, of which he was the first chair until 1987, and which he conceived and set up as a department including all the academic staff and the related teaching and research activities in computer and control of the University. He has been a visiting professor in universities and research centers in many different coun-

tries. From 1995 to 1998 he has been the President of the Italian Association of Researchers in Automatic Control, and from 1983 to present the coordinator of several national research projects in the field of automatic control. He is the author of many research papers spanning various diversified fields such as radar signal synthesis and processing, radio communications, system identification, adaptive, predictive, switching supervisory control, and detection of performance degradation in feedback-control systems. He is the author of a book Optimal, Predictive, and Adaptive Control (Prentice-Hall, 1995). He is an editor of the following journals: European Journal of Control; International Journal of Adaptive Control and Signal Processing, Wiley; and IEE Proceedings—Control Theory and Applications. Dr. Mosca is the Italian NMO representative in International Federation of Automatic Control (IFAC ). He has been a Council member of European Union Control Association (EUCA) until 1998, and from 1996 to 2002 a Council member of IFAC. In 2001, he was awarded “honoris causa” the Doctor in Information Engineering degree from the Universidade Tecnica de Lisboa, Lisbon, Portugal.