# Patent

## Title
Efficient digital filter design tool for approximating an FIR filter with a low-order linear-phase IIR filter

## Abstract

A method and apparatus for designing low-order linear-phase IIR filters is disclosed. Given an FIR filter, the method utilizes a new Krylov subspace projection method, called the rational Arnoldi method with adaptive orders, to synthesize an approximated IIR filter with small orders. The method is efficient in terms of computational complexity. The synthesized IIR filter can truly reflect essential dynamical features of the original FIR filter and indeed satisfies the design specifications. In particular, the linear-phase property is stilled remained in the passband.

## Background
BACKGROUND OF THE INVENTION
1. Field of the Invention
The present invention relates to a method of designing digital filters, and more particularly, to a method of designing low-order linear-phase IIR filters for approximating an FIR filter.
2. Description of Related Art
Digital filters for addressing digital signal processing have been applied widely in commercial electronic products, such as compact disk players, television sets and the like. In order to deal with real-time signals, the designs should be considered low computational costs to output signals efficiently. Also, since the distortion-free transmission of-the waveforms in the passband is very important in signal processing, the filters should contain a linear-phase characteristic, i.e., constant group delay.
There are two classes of digital filters. One is a finite-duration impulse response (FIR) filters and the other is an infinite-duration impulse response (IIR) filters. The main advantages of the FIR filters are exactly linear phase and guaranteed stable. However, when the specification being very rigorous, the resulting FIR filter is usually with higher orders, which may require more hardware components and lower the operational speed. Conversely, the IIR filters are useful for large-scale or high-speed designs but they do not have exactly linear phase and can not guarantee to be stable.
One way to synthesize the IIR filters with linear phase in passband is to solve the rational approximation problem directly. These methods, for example, include Pade approximation, linear programming, nonlinear programming, multiple criterion optimization and eigenfilter approach.
Another way is called the indirect approach. It will be composed of three steps, as shown in FIG. 1. Step 1 receives and stores the design specifications of a digital filter.
These specifications are required in the frequency-domain in terms of the desired magnitude and phase response of the filter. Then, a linear-phase FIR filter, which meets design specifications, will be designed in step 2. The order and the coefficients of the FIR filter can be obtained using the conventional methods such as the frequency-sampling design technique, the window design technique and the optimal equiripple design technique. Finally, a lower-order IIR filter will be obtained using filter approximation techniques in step 3. It should be mentioned that the special attentions shall be paid on this indirect approach. The resulting IIR filter must be ensured to capture the linear-phase response of the original FIR filter in the passband.
In recent years, several linear-phase IIR filter design techniques have been emerged for this purpose. Generally speaking, two distinct methods have be proposed: (1) Grammian-Based Methods: including the balanced truncation method and impulse response grammian method and (2) Optimal Approximation Methods: including the least-square approximations and the H2 norm approximation. Although satisfactory results have been reported, computational complexity of these methods are still quite expensive.
SUMMARY OF THE INVENTION
The main objective of the present invention is to provide a water-preventing grommet for a pull chain switch, which efficiently keeps water out of the pull chain switch to avoid malfunction.
The secondary objective of the present invention is to provide a water-preventing grommet for a pull chain switch, which smoothes operations of the pull chain switch.
To achieve the objectives, the method of approximating an FIR filter with low-order linear-phase IIR filters by the rational Arnoldi algorithm with adaptive orders in accordance with the present invention contains the following steps:

    
    
        (a) initialize the first vector of the Krylov sequence for each expansion point;
        (b) in the jth iteration of the algorithm, choose an expansion frequency such that the frequency gives the greatest difference between the (j+1)st-order output moment of the original FIR filter H(z) and that of the lower-order IIR filter Ĥ(z);
        (c) after the chosen expansion point in jth iteration being determined, the single-point Arnoldi method applied at the expansion point to generate the new orthnormal vector;
        (d) determine the new residual at each expansion point for the next iteration; wherein the resulting orthogonal projection matrix is output after giving the total iteration number of the algorithm.
    
    





BRIEF DESCRIPTION OF THE DRAWINGS
FIG. 1 illustrates an indirect approach to design a low-order linear-phase IIR filter;
FIG. 2 illustrates the design flow of the low-order linear-phase IIR filters;
FIG. 3 illustrates the typical design specifications of a low-pass filter;
FIG. 4 shows the detail flow of the rational Arnoldi method with adaptive orders; and
FIGS. 5A-7C show the bode plots of the magnitude, the error in magnitude, and the phase of the original FIR filters and the low-order IIR filters.



DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT
FIG. 2 shows the design flow of a low-order linear-phase IIR filter in the present invention in which includes the steps of receiving and storing the design specifications in step 1, designing an FIR filter satisfying the design specifications and saving the order and coefficients in step 2, establishing the state space matrices {A,b,c} in step 3, performing the rational Arnoldi method with adaptive orders and produce the orthogonal projection matrix V in step 4, and generating the corresponding low-order linear-phase IIR filter, which can approximate the original FIR filter and satisfy the design specifications in step 5.
FIG. 3 illustrates typical design specifications of a low-pass filter in step 1, where the band [0,ωp] (unit 34) is called the passband and δ1 (unit 36) is the acceptable tolerance (or ripple) in the passband, the band [ωs,π] (unit 38) is called the stopband and δ2 (unit 40) is the corresponding tolerance (or ripple), and the band [ωp,ωs] (unit 42) is called the transition band. Rp is the passband ripple in dB, where Rp=−20 log10 [(1−δ1)/(1+δ1)]. As is the stopband attenuation in dB, where As=−20 log10 [δ2/(1+δ1)]. Notably, either {δ1,δ2} or {Rp, As} is required to be stored in step 1.
Suppose that an FIR filter has been designed to satisfy the design specifications in step 2. Let H(z)=Σi=0nhiz−i be the causal FIR filter with length n+1. A state-space realization of H(z) in step 3 can be described as




  
    
      
        
          
            
              
                
                  x
                  ⁡
                  
                    (
                    
                      k
                      +
                      1
                    
                    )
                  
                
                =
                
                  
                    Ax
                    ⁡
                    
                      (
                      k
                      )
                    
                  
                  +
                  
                    bu
                    ⁡
                    
                      (
                      k
                      )
                    
                  
                
              
            
          
          
            
              
                
                  y
                  ⁡
                  
                    (
                    k
                    )
                  
                
                =
                
                  
                    
                      c
                      T
                    
                    ⁢
                    
                      
                    
                    ⁢
                    
                      x
                      ⁡
                      
                        (
                        k
                        )
                      
                    
                  
                  +
                  
                    
                      h
                      0
                    
                    ⁢
                    
                      u
                      ⁡
                      
                        (
                        k
                        )
                      
                    
                  
                
              
            
          
        
        ,
      
    
    
      
        (
        1
        )
      
    
  
  
    
      where
    
    
      
        
      
    
  
  
    
      
        
          
            
              
                A
                =
                
                  [
                  
                    
                      
                        0
                      
                      
                        0
                      
                      
                        ⋯
                      
                      
                        0
                      
                      
                        0
                      
                    
                    
                      
                        1
                      
                      
                        0
                      
                      
                        ⋯
                      
                      
                        0
                      
                      
                        0
                      
                    
                    
                      
                        0
                      
                      
                        1
                      
                      
                        ⋯
                      
                      
                        0
                      
                      
                        0
                      
                    
                    
                      
                        ⋮
                      
                      
                        ⋮
                      
                      
                        ⋮
                      
                      
                        ⋮
                      
                      
                        ⋮
                      
                    
                    
                      
                        0
                      
                      
                        0
                      
                      
                        ⋯
                      
                      
                        1
                      
                      
                        0
                      
                    
                  
                  ]
                
              
              ,
            
          
          
            
              
                b
                =
                
                  [
                  
                    
                      
                        1
                      
                    
                    
                      
                        0
                      
                    
                    
                      
                        ⋮
                      
                    
                    
                      
                        ⋮
                      
                    
                    
                      
                        0
                      
                    
                  
                  ]
                
              
              ,
            
          
          
            
              
                c
                =
                
                  [
                  
                    
                      
                        
                          h
                          1
                        
                      
                    
                    
                      
                        
                          h
                          2
                        
                      
                    
                    
                      
                        ⋮
                      
                    
                    
                      
                        ⋮
                      
                    
                    
                      
                        
                          h
                          n
                        
                      
                    
                  
                  ]
                
              
              ,
            
          
        
      
    
    
      
        (
        2
        )
      
    
  




and AεRn×n, bεRn, cεRn. The transfer function H(z) can also be expressed as H(z)=cTX(z)+h0=cT(zIn−A)−1b+h0. Our problem formulation is to find a lower-order IIR filter Ĥ(z), which satisfies the same specifications in step 1 as the original FIR filter H(z) and maintains a linear-phase response in the passband.

The way in the invention is to find an optimal IIR filter by using orthogonal projection of the original FIR filter. By matching some characteristics of the original FIR filter, the resulting orthonormal matrix V can be generated in step 4. The lower-order IIR filter Ĥ(z) can be constructed using the orthonormal projection x(k)=V{circumflex over (x)}(k). In such a situation, the parameters of the IIR filter can be defined by the following congruence transformation in step 5,

Â=VTAV, {circumflex over (b)}=VTb, and ĉ=VTc.  (3)

It can be shown that the matrix VT AV is always stable as long as (1) and matrix A is stable, (2) VTV=I. Thus, the stability of the lower-order IIR filter generated by Eq. (3) is guaranteed.

Pade Approximation and Moment Matching

The basis theory of the method in the invention is the multi-point Pade approximation, or so called the multi-point moment matching, to obtain a low-order IIR filter. Expanding X(z) in power series about various frequencies {z1, z2, . . . , zî}, where each zi=ejωiεC and 0≦ωi≦π, we have




  
    
      
        
          
            X
            ⁡
            
              (
              z
              )
            
          
          =
          
            
              ∑
              
                j
                =
                0
              
              ∞
            
            ⁢
            
              
                
                  X
                  
                    (
                    j
                    )
                  
                
                ⁡
                
                  (
                  
                    z
                    i
                  
                  )
                
              
              ⁢
              
                
                  (
                  
                    z
                    -
                    
                      z
                      i
                    
                  
                  )
                
                j
              
            
          
        
        ,
      
    
    
      
        (
        4
        )
      
    
  
  
    
      where
    
    
      
        
      
    
  
  
    
      
        
          
            
              X
              
                (
                j
                )
              
            
            ⁡
            
              (
              
                z
                i
              
              )
            
          
          =
          
            
              
                [
                
                  -
                  
                    
                      (
                      
                        
                          
                            z
                            i
                          
                          ⁢
                          
                            I
                            n
                          
                        
                        -
                        A
                      
                      )
                    
                    
                      -
                      1
                    
                  
                
                ]
              
              j
            
            ⁢
            
              
                (
                
                  
                    
                      z
                      i
                    
                    ⁢
                    
                      I
                      n
                    
                  
                  -
                  A
                
                )
              
              
                -
                1
              
            
            ⁢
            b
          
        
        ,
      
    
    
      
        (
        5
        )
      
    
  
  
    
      
        
          
            
              
                
                  H
                  
                    (
                    j
                    )
                  
                
                ⁡
                
                  (
                  
                    z
                    i
                  
                  )
                
              
              =
              
                
                  c
                  T
                
                ⁢
                
                  
                
                ⁢
                
                  
                    X
                    
                      (
                      j
                      )
                    
                  
                  ⁡
                  
                    (
                    
                      z
                      i
                    
                    )
                  
                
              
            
          
          
            
              
            
          
          
            
              
                (
                
                  j
                  >
                  0
                
                )
              
              ,
            
          
        
      
    
    
      
        
      
    
  
  
    
      
        
          
            H
            
              (
              0
              )
            
          
          ⁡
          
            (
            
              z
              i
            
            )
          
        
        =
        
          
            
              c
              T
            
            ⁢
            
              
            
            ⁢
            
              
                X
                
                  (
                  0
                  )
                
              
              ⁡
              
                (
                
                  z
                  i
                
                )
              
            
          
          +
          
            
              h
              0
            
            .
          
        
      
    
    
      
        
      
    
  




X(j)(zi) is called the jth-order system moment of X(z); H(j)(zi) represents the jth-order output moment of H(z) at zi. Notably, if î=1, Eq. (4) is indeed the conventional Pade approximation. The objective is to find a q-order (q<n) IIR filter Ĥ(z)=ĉT(zIq−Â)−1{circumflex over (b)}+h0 such that H(j)(zi)=Ĥ(j)(zi) for j=0,1, . . . , ĵi−1 and i=1,2, . . . , î, where q=Σi=1îĵi.

It shall be mentioned that moment calculations can be obtained analytically by exploring special characteristics of matrices A and b in eq. (2). For each zi, (ziIn−A)−1b and (ziIn−A)−1 can be derived analytically as the following formulas:




  
    
      
        
          
            
              (
              
                
                  
                    z
                    i
                  
                  ⁢
                  
                    I
                    n
                  
                
                -
                A
              
              )
            
            
              -
              1
            
          
          =
          
            [
            
              
                
                  
                    1
                    /
                    
                      z
                      i
                    
                  
                
                
                  0
                
                
                  ⋯
                
                
                  0
                
                
                  0
                
              
              
                
                  
                    1
                    /
                    
                      z
                      i
                      2
                    
                  
                
                
                  
                    1
                    /
                    
                      z
                      i
                    
                  
                
                
                  ⋯
                
                
                  0
                
                
                  0
                
              
              
                
                  
                    1
                    /
                    
                      z
                      i
                      3
                    
                  
                
                
                  
                    1
                    /
                    
                      z
                      i
                      2
                    
                  
                
                
                  ⋯
                
                
                  0
                
                
                  0
                
              
              
                
                  ⋮
                
                
                  ⋮
                
                
                  ⋮
                
                
                  ⋮
                
                
                  ⋮
                
              
              
                
                  
                    1
                    /
                    
                      z
                      i
                      
                        n
                        -
                        1
                      
                    
                  
                
                
                  
                    1
                    /
                    
                      z
                      i
                      
                        n
                        -
                        2
                      
                    
                  
                
                
                  ⋯
                
                
                  
                    1
                    /
                    
                      z
                      i
                    
                  
                
                
                  0
                
              
              
                
                  
                    1
                    /
                    
                      z
                      i
                      n
                    
                  
                
                
                  
                    1
                    /
                    
                      z
                      i
                      
                        n
                        -
                        1
                      
                    
                  
                
                
                  ⋯
                
                
                  
                    1
                    /
                    
                      z
                      i
                      2
                    
                  
                
                
                  
                    1
                    /
                    
                      z
                      i
                    
                  
                
              
            
            ]
          
        
        ,
      
    
  
  
    
      
        
          
            
              (
              
                
                  
                    z
                    i
                  
                  ⁢
                  
                    I
                    n
                  
                
                -
                A
              
              )
            
            
              -
              1
            
          
          ⁢
          b
        
        =
        
          
            
              [
              
                
                  
                    
                      1
                      /
                      
                        z
                        i
                      
                    
                  
                  
                    
                      1
                      /
                      
                        z
                        i
                        2
                      
                    
                  
                  
                    ⋯
                  
                  
                    
                      1
                      /
                      
                        z
                        i
                        n
                      
                    
                  
                
              
              ]
            
            T
          
          .
        
      
    
  




Krylov Subspace and the Arnoldi Method

Explicitly computing moments usually yields numerically ill-conditioned problems. We adapt recent works about the Krylov space method to solve these problems. Given a square matrix ΨεCn×n and a vector ξεCn, the qth Krylov sequence

Kq(Ψ, ξ)≡span(ξ,Ψξ,Ψ2ξ, . . . ,Ψq−1ξ)

is a sequence of q column vectors and the corresponding column space is called the qth Krylov subspace. Set Ψ=(ziIn−A)−1 and ξ=(ziIn−A)−1b. It has been shown that the Krylov subspace Kq(Ψ,ξ) is indeed spanned by the system moments X(j)(zi) for j=0,1, . . . , q−1. The Arnoldi method, a kind of Krylov subspace methods, is employed to generate an orthonormal matrix Vq that spans the same subspace as the Krylov subspace Kq(Ψ,ξ). As a result, the guaranteed stable IIR filter can be constructed by substituting Vq into Eq. (3).

The Arnoldi method arises from the Hessenberg reduction A=VHVT for eigenvalue calculations. It has the advantage that it can be terminated part-way and leaving one with a partial reduction to a Hessenberg form. The process is exploited to form iterative algorithms. During the iteration process, an upper Hessenberg matrix HqεCq×q is generated that satisfies the following relationship:

ΨVq=VqHq+hq+1,qvq+1eqT and v1=ξ/∥ξ∥,  (6)

where eq is the qth unit vector in Rq. The vector vq+1 satisfies a (q+1)-term recurrence relation, involving itself and the preceding Krylov vectors. A new orthonormal vector vq+1 can be generated using the modified Gram-Schmidt orthogonalization technique.

The Rational Arnoldi Method

Generally speaking, the accuracy of the Pade approximation based methods is lost away from the expansion point more rapidly as the eigenvalues of the FIR filter approach the expansion frequency. A rational Arnoldi (RA) method, which uses multiple expansion points, was developed to overcome this difficulty. The straightforward way for multi-point moment matching applications is to apply the Krylov subspace algorithm at various expansion frequencies. This is the so-called rational Krylov algorithm. Basically, this algorithm is a generalization of the shifted-and-inverted Arnoldi algorithm. To simplify the developments, the number of the matched moments of the lower-order IIR filter at each expansion point is assumed to be fixed. Formally, let Z={z1, z2, . . . , zî} represent the set of predetermined expansion frequencies. Let J={ĵ1, ĵ2, . . . , ĵî} be the set of the number of the matched moments at each corresponding frequency. The rational Arnoldi method will generate a lower-order IIR filter Ĥ(z), which matches q-order (q=Σi=1îĵi) moments of the FIR filter, H(z), at the expansion points zi, i=1,2, . . . , î.
Implementing the rational Arnoldi method is equivalent to implement the Arnoldi method ĵi times at î expansion frequencies. That is, the first ĵ1 iterations correspond to the expansion frequency z1 and the next ĵ2 iterations are associated with z2, and so on. Each Arnoldi iteration generates ĵi orthonormal vectors. Then, Vq=└v1 v2 . . . vq┘ is the desired orthonormal matrix generated from a union Krylov space at various expansion points, as stated by

Kq=span(X(0)(z1), . . . ,X(ĵ1−1)(z1), . . . ,X(0)(zî), . . . ,X(ĵî−1)(zî)).

Once the orthonormal matrix Vq has been formed by applying the rational Arnoldi method and the lower-order IIR filter can be obtained using the congruence transformation.
The Rational Arnoldi Method with Adaptive Orders
Selecting a set of expansion points zi for i=1,2, . . . , î and the number of matched moments ĵi about each zi is by no means trivial. For simplicity, the expansion points zi for i=1,2, . . . , î are determined in advance using engineering heuristics or experimental measurements over a specified frequency range. This invention describes an intelligent scheme for choosing multiple expansion points in each of the iterations.
Suppose that H(j)(zi)=Ĥ(j)(zi) for j=0,1, . . . , ĵi−1 and i=1,2, . . . , î after q iterations of the rational Arnoldi algorithm. However, the ĵith-order output moments H(ĵi)(zi)=Ĥ(ĵi)(zi) can not be guaranteed. The concept that underlies the rational Arnodli method with adaptive orders is to select an expansion point zi*q+1 in the (q+1)st iteration. Hence, the resulting (q+1)st-order IIR filter yields the greatest moment improvement |H(ĵi)(zi)−Ĥ(ĵi)(zi)| of the qth-order IIR filter as zi=zi*q+1. The moment errors can be directly obtained in the new iteration without explicitly calculating system moments.
FIG. 4 displays the detail flow of the rational Arnoldi method with adaptive orders in step 4 in FIG. 2.
Step 1, in FIG. 4, initializes the first vector k(0)(zi)=(ziIn−A)−1b of the Krylov sequence for each expansion point zi, where iε{1, . . . , î}. Since the lower-order IIR filter and the orthonormal matrix are not yet determined, the residue r(0)(zi) for each zi is set to k(0)(zi). The normalization coefficient about each zi, hπ(zi), is initialized to be one. Step 2, in FIG. 4, begins the iterations and sets j=1.
Step 3, in FIG. 4, chooses an expansion frequency zi such that zi gives the greatest difference between the (j+1)st-order output moment of the original FIR filter H(z) and that of the lower-order IIR filter Ĥ(z), that is,

maxziεZ|H(j+1)(zi)−Ĥ(j+1)(zi)|=maxziεZ|hπ(zi)cTr(j−1)(zi)|.

Ĥ(j+1)(zi) is the (j+1)st-order output moment of the lower-order IIR filter Ĥ(z), which is yielded using the congruence transformation matrix Vj−1 (j>1) and matches j-order output moments of H(z) at zi. The chosen expansion frequency in the jth iteration is called zi*j.

After choosing the expansion point zi*j in the determined jth iteration, the single-point Arnoldi method is applied at the expansion point zi*j (unit 52), which contains steps 4 and 5, as shown in FIG. 4. Step 4, in FIG. 4, generates the new orthnormal vector vj and the vector is incorporated into the orthnormal matrix Vj−1. The normalization coefficient hπ(zi)=Πj∥r(j−1)(zi)∥ when zi is selected in the jth iteration.
Step 5, in FIG. 4, determines the new residual r(j)(zi) at each expansion point zi. The calculation involves a projection with the new orthonormal matrix Vj. The next vector k(j)(zi*q+1) at the frequency zi*q+1 must be updated to enable further matching of the output moment in the (j+1)st iteration. Since no improvement is obtained at the other unselected frequency zi, the vector k(j)(zi) at frequency zi in the current iteration remains k(j−1)(zi), which was obtained in the preceding iteration. Reset j=j+1 in step 6 and judge if j≦q in step 7, as shown in FIG. 4. Finally, the resulting orthogonal projection matrix Vq is generated in step 8 in FIG. 4.
The resulting orthnormal matrix Vq should be real to ensure that real system matrices of the lower-order IIR filter are generated if the complex expansion frequencies are used. First, all column vectors in Vq are divided into the real part Vr and the imaginary part Vi. Second, a reduced QR factorization of [Vr Vi] is performed to yield a new orthogonal matrix Vq. The moment matching property of the resulting lower-order IIR filter by the new and real Vq is also preserved.
The details of the algorithm are outlined as follows. The vector Z includes î expansion points, q is the total number of iterations and Vq is the resulting orthonormal matrix.










Adaptive Rational Arnoldi (input: A,b,c,Z,q; output: Vq)

















(1): /* Initialize */



1 for each zi ∈ Z do



2   k(0)(zi) := (ziIn−A)−1b, r(0)(zi) := k(0)(zi)



3   hπ(zi) := 1



4 end for



(2): /* Begin the Iterations */



5 for j = 1, 2, . . ., q do



    (2.1) /* Select the Expansion Frequency with the Maximum



       Output Moment Error*/



6   Choose zi ∈ Z as the i giving maxi|hπ(zi)cTr(j−1)(zi)|



7   set zi*j be the expansion frequency in the jth iteration



    (2.2) /*Generate the Orthonormal Vector at zi*j */



8   hj,j−1(zi*j) := ||r(j−1)(zi*j)||



9   vj = r(j−1)(zi*j)/hj,j−1(zi*j)



10   hπ(zi*j) := hπ(zi*j)·hj,j−1(zi*j)



    (2.3) /* Update the Residue r(j)(zi) for the Next Iteration */



11    for each zi ∈ Z do



12     if (zi == zi*j) then k(j)(zi*j) := −(ziIn−A)−1vj



13     else k(j)(zi) := k(j−1)(zi)



14     end if



15     r(j)(zi) := k(j)(zi)



16       for t = 1, 2, . . ., j do



17         ht,j(zi) := vtHr(j)(zi)



18         r(j)(zi) := r(j)(zi)−ht,j(zi)vt



19       end for



20    end for



21 end for



22 Vq = [V1 V2 . . . Vq]










Some properties of the method of approximating an FIR filter by low-order IIR filters in the invention are summarized as follows.
(1) Exact expression of output moment errors: suppose that the output moments of the original FIR filter and those of the lower-order IIR filter are matched, that is, H(j)(zi)=Ĥ(j)(zi) for j=0,1, . . . , ĵi−1 and i=1,2, . . . , î. The system matrices of the lower-order IIR filter are generated by the congruence transformation with the orthonormal matrix Vq using the algorithm, where q=Σi=1îĵi. The magnitude error between the ĵith-order moments H(ĵi)(zi) and Ĥ(ĵi)(zi) at each expansion point zi can be expressed as follows:

|H(ĵi)(zi)−Ĥ(ĵi)(zi)|=|hπcTr(ĵi−1)(zi)|,  (10)

where hπ(zi)=Πj∥r(j−1)(zi)∥.

(2) Moment matching can still be preserved.
(3) In the first iteration in the rational Arnodli algorithm with adaptive orders, step (2.2) is to choose ziεZ such that max(|cT(ziIn−A)−1b|)=max(|H(zi)|). This is equivalent to find out the expansion frequency with the maximum magnitude in the output frequency response.
(4) Implementation issues of digital filters: the present invention also provides several heuristics of selecting expansion frequencies in advance for the proposed rational Arnoldi method. Generally speaking, the complex expansion points {z1, z2, . . . , zî} will be recommended, where each zi=ejωiεC and 0≦ωi≦π. Then the frequency responses of the lower-order IIR filters at these points can be the same as those of the original FIR filter. Nevertheless, if real expansion points can be selected, the computational complexity of yielding approximate IIR filters can be further reduced. The following guidelines are provided:
(a) Low-pass/high-pass filters: the proposed method with the expansion point ω1=0 performs well over the low frequency range of responses. For high-pass filter designs, the special structures of state-space matrices may be used to present the duality between low-pass and high-pass filters. Let Ā=−A, b=b, c=c, and h0=−h0,




  
    
      H
      _
    
    ⁡
    
      (
      z
      )
    
  
  =
  
    
      
        
          
            
              
                c
                _
              
              T
            
            ⁡
            
              (
              
                
                  zI
                  n
                
                -
                
                  A
                  _
                
              
              )
            
          
          
            -
            1
          
        
        ⁢
        
          b
          _
        
      
      +
      
        
          h
          _
        
        0
      
    
    =
    
      
        ∑
        
          i
          =
          0
        
        n
      
      ⁢
      
        
          
            (
            
              -
              1
            
            )
          
          
            i
            +
            1
          
        
        ⁢
        
          h
          i
        
        ⁢
        
          
            z
            
              -
              i
            
          
          .
        
      
    
  




If H(z) presents a high-pass filter, then H(z) will be a low-pass filter, and a vice versa. Likewise, the expansion point ω1=0 is chosen to perform the Arnoldi algorithm. If the corresponding orthonormal matrix Vq is obtained, then the high-pass IIR filter, which satisfies the same specifications as the original FIR filter, can be constructed as follows:

Â= VqTA Vq, {circumflex over (b)}= VqTb, and ĉ= VqTc.

(b) Band-pass/band-stop filters: experimental results indicate that the passband edge and stopband edge frequencies are appropriate candidate expansion points in meeting the specifications of the design. Other expansion points with uniform spacing are also recommend to be selected.
DESIGN EXAMPLES
Three example filters are used to justify the proposed approach. Table 1 describes specifications of a low-pass filter, a high-pass filter, and a band-pass filter. The command remez in Matlab was used to design the FIR filters by the optimal equiripple technique. Table 2 lists the corresponding orders. Then, the approximate low-order IIR filters were generated by the proposed method and the balanced realization method (BAL). Table 2 shows the reduced orders and the expansion points used by the two methods. FIGS. 5A-7C display the bode plots of the magnitude, the error in magnitude, and the phase of the original FIR filters and the low-order IIR filters. In FIGS. 5A-7C, the responses of the original FIR filters are represented as thin solid lines. Those of the IIR filters, determined by the proposed method, are represented as thick solid lines -, and those determined by BAL method are plotted as thick dashed lines - -. The responses in the passband of the IIR filters are indistinguishable from those of the original FIR filters, independently of which the model reduction method is used. Simulation results imply that the performance of the proposed method is similar to that of the BAL method in the passband. The resulting lower-order IIR filters can actually preserve the linear-phase response of the original FIR filters. Nevertheless, in terms of computational efficiency, the Kylov subspace based methods generally outperform the BAL method.







TABLE 1







Filter design specifications









Specifications











Low-Pass
High-Pass
Band-Pass














Maximum passband attenuation
3
1
1


(dB)


Maximum stopband attenuation
40
35
35


(dB)


Lower passband edge (rad/s)
0
0.85π
0.3255π


Upper passband edge (rad/s)
0.285π
  1π
0.3755π


Lower stopband edge (rad/s)
0.353π
0
0.6655π


Upper stopband edge (rad/s)
   1π
0.78π
0.7155π
















TABLE 2







Matched orders and expansion points for IIR filter designs











Low-Pass
High-Pass
Band-Pass















FIR
Order
41
41
58


IIR
Order
24
17
36











Expansion
ω = 0
ω = {0.4π, 0.6π}



points









CONCLUSIONS
A rational Arnoldi method with adaptive orders for approximating FIR filters by low-order linear-phase IIR filters has been proposed. The developed method is very efficient in terms of computational complexity. Meanwhile, the lower-order IIR filter can truly reflect the dynamical features of the FIR filter and satisfies the original design specifications.
Although the invention has been explained in relation to its preferred embodiment, it is to be understood that many other possible modifications and variations can be made without departing from the spirit and scope of the invention as hereinafter claimed.

## Claims

1. A method of approximating an FIR filter with low-order linear-phase IIR filters by the rational Arnoldi algorithm with adaptive orders containing the following steps:
a) initialize the first vector of the Krylov sequence for each expansion point;
b) in the jth iteration of the algorithm, choosing an expansion frequency wherein the heuristics of selecting expansion frequencies in advance for the proposed rational Arnoldi method we given by
(a) low-pass filters: the proposed method with the expansion point ω1=0;
(b) high-pass filters: the special structures of state-space matrices used to present the duality between low-pass and light pass filters; let state matrices become Ā=−A, b=b, c=c, and h0=−h0, the expansion point ω1=0 chosen to perform the Arnoldi algorithm; when the corresponding orthonormal matrix Vq is obtained and then the high-pass IIR filter, which satisfies the same specifications as the original FIR filter; and
(c) band-pass/band-stop filters: the passband edge and the stopband edge frequencies being the appropriate candidate expansion points in meeting the specifications of the design, and other expansion points with uniform spacing recommended to be selected


such that the frequency gives the greatest difference between the (j+1)st-order output moment of the original FIR filter H(z) and that of the lower-order IIR filter Ĥ(z) wherein the expression of output moment errors between the ĵith-order moments H(ĵi)(zi) and Ĥ(ĵi)(zi) at each expansion point zi are expressed as follows:



  
    
      
        
          
          
            
              
                H
                
                  (
                  
                    
                      j
                      ^
                    
                    i
                  
                  )
                
              
              ⁡
              
                (
                
                  z
                  i
                
                )
              
            
            -
            
              
                
                  H
                  ^
                
                
                  (
                  
                    
                      j
                      ^
                    
                    i
                  
                  )
                
              
              ⁡
              
                (
                
                  z
                  i
                
                )
              
            
          
          
        
        =
        
          
          
            
              h
              π
            
            ⁢
            
              c
              T
            
            ⁢
            
              
                r
                
                  (
                  
                    
                      
                        j
                        ^
                      
                      i
                    
                    -
                    1
                  
                  )
                
              
              ⁡
              
                (
                
                  z
                  i
                
                )
              
            
          
          
        
      
      ,
      
        

      
      ⁢
      where
    
    ⁢
    
      
    
  




  
    
      
    
    ⁢
    
      
        
          h
          π
        
        ⁡
        
          (
          
            z
            i
          
          )
        
      
      =
      
        
          ∏
          j
        
        ⁢
        
          
        
        ⁢
        
          
          
            
            
              
                r
                
                  (
                  
                    j
                    -
                    1
                  
                  )
                
              
              ⁡
              
                (
                
                  z
                  i
                
                )
              
            
            
          
          
        
      
    
  



is the normalization coefficient when an expansion frequency zi is selected in the jth iteration; vector c contains the last n impulse response coefficients of a FIR filter with length n+1; and r(j−1)(zi) is the residual vector in the (j−1)st iteration of the disclosed adaptive rational Arnoldi algorithm at the expansion frequency zi;
c) after the choosing the expansion point in jth iteration being determined, the single-point Arnoldi method applied at the expansion point to generate the new orthnormal vector; and
d) determine a new residual at each expansion point for next iteration;

whereby, after the giving total iteration number of the algorithm, outputting the resulting orthogonal projection matrix.