Skip to content Skip to navigation

Connexions

You are here: Home » Content » Direct Frequency Domain IIR Filter Design Methods

Navigation

Content Actions
  • Download module PDF
  • Add to ...
    Add the module to:
    • My Favorites
    • A lens
    • An external social bookmarking service
    • My Favorites (What is 'My Favorites'?)
      'My Favorites' is a special kind of lens which you can use to bookmark modules and collections directly in Connexions. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need a Connexions account to use 'My Favorites'.
    • A lens (What is a lens?)

      Definition of a lens

      Lenses

      A lens is a custom view of Connexions content. You can think of it as a fancy kind of list that will let you see Connexions through the eyes of organizations and people you trust.

      What is in a lens?

      Lens makers point to Connexions materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

      Who can create a lens?

      Any individual Connexions member, a community, or a respected organization.

    • External bookmarks
  • E-mail the author

Recently Viewed

Direct Frequency Domain IIR Filter Design Methods

Module by: C. Sidney Burrus

The preceding design methods have been based on designing an analog prototype filter and then converting it to a digital filter. This approach is appropriate for the class of approximations where analytical solutions are possible, but not for many others. In the remaining part of this chapter, methods will be developed that directly design the desired digital filter. Most approaches are extensions of methods used for FIR filters, but they are more complicated for the IIR case where rational approximation is being performed rather than polynomial approximation.

In this section a frequency-sampling design method is developed such that the frequency response of the IIR filter will pass through the given samples of a desired response. Since an IIR filter cannot have linear phase, the sampled response must contain both magnitude and phase. The extension of the frequency- sampling method to a LS-error approximation is not as simple as for the FIR filter. The method presented in this section uses a criterion based on the equation error rather than the more common error between the actual and desired frequency responseEntry 1. Nevertheless, it is a useful noniterative design method. Finally, a general discussion of iterative design methods for LS-frequency response error is given.

Frequency-Sampling Design of IIR Filters

The method for calculating samples of the frequency response of an IIR filter presented in the section on Properties of IIR Filters can be reversed to design a filter much the same way it was for the FIR filter using frequency sampling. The z-transform transfer function for an IIR filter is given by

H ( z ) = B ( z ) A ( z ) = b 0 + b 1 z - 1 + ⋯ + b M z - M 1 + a 1 z - 1 + ⋯ + a N z - N . H ( z ) = B ( z ) A ( z ) = b 0 + b 1 z - 1 + ⋯ + b M z - M 1 + a 1 z - 1 + ⋯ + a N z - N . (1)

The frequency response of the filter is given by setting z=e-jωz=e-jω. Using the notation

H ( ω ) = H ( z ) | z = e - j ω . H ( ω ) = H ( z ) | z = e - j ω . (2)

Equally-spaced samples of the frequency response are chosen so that the number of samples is equal to the number of unknown coefficients in (Equation 1). These L+1L+1 = M+N+1M+N+1 samples of this frequency response are given by

H k = H ( ω k ) = H ( 2 π k L + 1 ) H k = H ( ω k ) = H ( 2 π k L + 1 ) (3)

and can be calculated from the length-(L+1)(L+1)) DFTs of the numerator and denominator.

H k = DFT { b n } DFT { a n } = B k A k H k = DFT { b n } DFT { a n } = B k A k (4)

where the indicated division is term-by-term division for each value of kk. Multiplication of both sides of (Equation 4) by AkAk gives

B k = H k A k B k = H k A k (5)

If the length-(L+1)(L+1) inverse DFT of HkHk is denoted by the length- (L+1)(L+1) sequence hnhn, equation (Equation 5) becomes cyclic convolution which can be expressed in matrix form by

b 0 b 1 ⋮ b M 0 ⋮ 0 = h 0 h L h L - 1 ⋯ h 1 h 1 h 0 h L h 2 h 1 h 0 ⋮ ⋮ h L ⋯ h 0 1 a 1 ⋮ a N 0 ⋮ 0 b 0 b 1 ⋮ b M 0 ⋮ 0 = h 0 h L h L - 1 ⋯ h 1 h 1 h 0 h L h 2 h 1 h 0 ⋮ ⋮ h L ⋯ h 0 1 a 1 ⋮ a N 0 ⋮ 0 (6)

Note that the hnhn in (Equation 6) are not the impulse response values of the filter as used in the FIR case. A more compact matrix notation of is

b 0 = H a 0 b 0 = H a 0 (7)

where HH is (L+1)(L+1) by (L+1)(L+1), bb is length-(M+1)(M+1), and aa is length-(N+1)(N+1). Because the lower L-NL-N terms of the right-hand vector of (Equation 6) are zero, the HH matrix can be reduced by deleting the right-most L-NL-N columns to give H0H0 which causes (Equation 7) to become

b 0 = H 0 a b 0 = H 0 a (8)

Because the first element of aa is unity, it is partitioned to remove the unity term and the remaining length-NN vector is denoted a*a*. The simultaneous equations represented by (Equation 8) are uncoupled by further partitioning of the HH matrix as shown in

b 0 = H 1 h 1 H 2 1 a * b 0 = H 1 h 1 H 2 1 a * (9)

where H1H1 is (M+1)(M+1) by (N+1)(N+1), h1h1 is length-(L-M)(L-M), and H2H2 is (L-M)(L-M) by NN. The lower (L-M)(L-M) equations are written

0 = h 1 + H 2 a * 0 = h 1 + H 2 a * (10)

or

h 1 = - H 2 a * h 1 = - H 2 a * (11)

which must be solved for a*a*. The upper M+1M+1 equations of (10) are written

b = H 1 a b = H 1 a (12)

which allows the calculation of bb.

If L=N+ML=N+M, H2H2 is square. If H2H2 is nonsingular, (Equation 11) can be solved exactly for the denominator coefficients in a*a*, which are augmented by the unity term to give aa. From (Equation 12), the numerator coefficients in bb are found. If H2H2 is singular Entry 16, (Reference) and there are multiple solutions, a lower order problem can be posed. If there are no solutions, the approximation methods must be used.

Note that any order numerator and denominator can be prescribed. If the filter is in fact an FIR filter, aa is unity and a*a* does not exist. Under these conditions, (Equation 12) states that bn=hnbn=hn, which is one of the cases of FIR frequency sampling covered Entry 19. Also note that there is no control over the stability of the filter designed by this method.

ummary

In this section, an interpolation design method was developed and analyzed. Use of the DFT converted the frequency- domain specifications to the time domain. A matrix partitioning allowed uncoupling the solution for the numerator from the solution of the denominator coefficients. The use of the DFT prevents the possibility of unequally spaced frequency samples as was possible for FIR filter design. The solution of simultaneous equations would allow unequal spacing which is not as troublesome as with the FIR filter because IIR filters are usually of lower order.

The frequency-sampling design of IIR filters is somewhat more complicated than for FIR filters because of the requirement that H2H2 be nonsingular. As for the FIR filter, the samples of the desired frequency response must satisfy the conditions to insure that hnhn are real. The power of this method is its ability to interpolate arbitrary magnitude and phase specification. In contrast to most direct IIR design methods, this method does not require any iterative optimization with the accompanying convergence problems.

As with the FIR version, because this design approach is an interpolation method rather than an approximation method, the results may be poor between the interpolation points. This usually happens when the desired frequency-response samples are not consistent with what an IIR filter can achieve. One solution to this problem is the same as for the FIR case Entry 19, the use of more frequency samples than the number of filter coefficients and the definition of an approximation error function that can be minimized. There is no simple restriction that will guarantee stable filters. If the frequency-response samples are consistent with an unstable filter, that is what will be designed.

Discrete Least-Squared Equation-Error IIR Filter Design

In order to obtain better practical filter designs, the interpolation scheme of the previous section is extended to give an approximation design method Entry 19. It should be noted at the outset that the method developed in this section minimizes an equation-error measure and not the usual frequency-response error measure.

The number of frequency samples specified, L+1L+1, will be made larger than the number of filter coefficients, M+N+1M+N+1. This means that H2H2 is rectangular and, therefore, (Equation 11) cannot in general be satisfied. To formulate an approximation problem, a length-(L+1)(L+1) error vector εε is introduced in (Equation 8) and (Equation 9) to give

b 0 = H 0 a + ε b 0 = H 0 a + ε (13)

Equation (Equation 11) becomes

h 1 - ε = - H 2 a * h 1 - ε = - H 2 a * (14)

where now H2H2 is rectangular with (L-M)>N(L-M)>N. Using the same methods as used to derive (Equation 11), the error εε is minimized in a least-squared error sense by the solution of the normal equations Entry 16

H 2 T h 1 = - H 2 T H 2 a * H 2 T h 1 = - H 2 T H 2 a * (15)

If the equations are not singular, the solution is

a * = - [ H 2 T H 2 ] - 1 H 2 T h 1 . a * = - [ H 2 T H 2 ] - 1 H 2 T h 1 . (16)

If the normal equations are singular, the pseudo-inverse Entry 16, (Reference) can be used to obtain a minimum norm or reduced order solution.

The numerator coefficients are found by the same techniques as before in (Equation 12)

b = H 1 a b = H 1 a (17)

which results in the upper M+1M+1 terms in εε being zero and the total squared equation error being minimum.

As is true for LS-error design of FIR filters, (Equation 15) is often numerically ill-conditioned and (Equation 16) should not be used to solve for a*a*. Special algorithms such as those used by Matlab and LINPACK Entry 18, Entry 7 should be employed.

The error εε defined in (Equation 13) can better be understood by considering the frequency-domain formulation. Taking the DFT of (Equation 13) gives

B k = H k A k + ε B k = H k A k + ε (18)

where εε is the error in trying to satisfy (Equation 8) when the equations are over-specified. This can be reformulated in terms of EE, the difference between the frequency response samples of the designed filter and the desired response samples, by dividing (Equation 8) by AkAk to give

E k = B k A k - H k = ε k A k E k = B k A k - H k = ε k A k (19)

EE is the error in the solution of the approximation problem, and εε is the error in the equations defining the problem. The usual statement of a frequency-domain approximation problem is in terms of minimizing some measure of EE, but that results in solving nonlinear equations. The design procedure developed in this section minimizes the squared error εε, thus only requiring the solution of linear equations. There is an important relation between these problems. Equation (Equation 19) shows that minimizing εε is the same as minimizing EE weighted by AA. However, AA is unknown until after the problem is solved.

Although this is posed as a frequency-domain design method, the method of solution for both the interpolation problem and the LS equation-error problem is the same as the time-domain Prony's method, discussed in Section 7.5 of reference Entry 19.

Numerous modifications and extensions can be made to this method. If the desired frequency response is close to what can be achieved by an IIR filter, this method will give a design approximately the same as that of a true least-squared solution-error method. It can be shown that ε=0↔E=0ε=0↔E=0. In some cases, improved results can be obtained by estimating AkAk and using that as a weight on εε to approximate minimizing EE. There are iterative methods based on solving (Equation 16) and (Equation 17) to obtain values for AkAk. These values are used as weights on εε to solve for a new set of AkAk used as a new set of weights to solve again for AkAkEntry 19Entry 26. We found this approach to converge slowly, but a recent paper using the log-magnitude Entry 14 was more successful. Other approaches are given in Entry 24, Entry 23, Entry 11. The solution of (Equation 16) and (Equation 17) is sometimes used to obtain starting values for other iterative optimization algorithms that need good starting values for convergence.

An interesting iterative design algorithm that can design to approximate complex or magnitude frequency responses has be recently proposed by Jackson Entry 11. A different approach to the same problem was posed by Soewito Entry 26, Entry 28.

To illustrate this design method a sixth-order lowpass filter was designed with 41 frequency samples to approximate. The magnitude of those less than 0.2 Hz is one and of those greater than 0.2 is zero. The phase was experimentally adjusted to result in a good magnitude response. The design was performed with Program 9 in the appendix of Entry 19 and the frequency response is shown in Figure 7-33 of Entry 19. Matlab programs have recently been written which are smaller and easier to understand than those in FORTRAN.

ummary

In this section an LS-error approximation method was posed to design IIR filters. By using an equation-error rather than a solution-error criterion, a problem resulted that required only the solution of simultaneous linear equations.

Like the FIR filter version, the IIR frequency sampling design method and the LS equation-error extension can be used for complex approximation and, therefore, can design with both magnitude and phase specifications.

If the desired frequency-response samples are close to what an IIR filter of the specified order can achieve, this method will produce a filter very close to what a true least-squared error method would. However, when the specifications are not consistent with what can be achieved and the approximating error is large, the results can be very poor and in some cases, unstable. It is particularly difficult to set realistic phase response specifications. With this method, it is even more important to have a design environment that will allow easy trial-and-error procedure.

Newly published works which will be discussed here are Entry 14, Entry 13, Entry 22, Entry 10, Entry 20, Entry 17, Entry 12. Other references can be found in Entry 19, Entry 14, Entry 6, Entry 26, Entry 28, Entry 5. The Matlab command invfreqz() which is an inverse to the freqz() command gives a similar or, perhaps, the same result as the method described in this note but uses a different formulation Entry 15, Entry 25.

more

Practical problems occur in the design of a filter to separate signals according to their energy. Because the energy content of a signal is the integral or sum of the square of the signal, a mean-squared-error measure is natural. Unfortunately, for the IIR filter design problem, the optimization procedure is nonlinear. This was pointed out in the last section where the equation error was used in order to have a linear problem.

Because of the nonlinear nature of the least-squared-error minimization, the method of solution becomes dependent on the desired frequency response, and therefore, there is no single method for design. The mean-squared error for magnitude approximation is defined as

q ( x ) = ∑ i = 0 L | H ( ω i ) | - | H d ( ω i ) | 2 q ( x ) = ∑ i = 0 L | H ( ω i ) | - | H d ( ω i ) | 2 (20)

where xx is a vector of filter parameters chosen to minimize q, and the error is sampled at L+1L+1 frequencies ωiωi. Steiglitz Entry 19 chose the parameter vector x to be the coefficients of a cascade structure in order to best fit an iterative optimization scheme. He applied a standard optimization algorithm, the Fletcher-Powell method, to the minimization of (Equation 19). Other methods which are more directly related to a squared-error measure can also be used.

Practical difficulties exist in solving this approximation problem. In some cases, local minima are found rather than the global minimum. In other cases, convergence of the minimization algorithm is slow or does not occur at all. Numerical problems can result from ill-conditioned equations, and there is no guarantee that the designed filter will be stable.

An important factor is the choice of a desired frequency- response function Hd(ω)Hd(ω) that does not result in the optimum approximation having a large error. This often means not having an abrupt discontinuity between the passband and stopband.

Another factor is the starting of the iterative optimization algorithm with a set of coefficients in xx that is close to the optimum. This can be accomplished by using the frequency sampling method to give a design that can be used to start a least-squares algorithm. Because the error defined in (Equation 19) is in terms of magnitudes, an unstable design can be converted to a stable one by moving the unstable pole at a radius of rr in the zz-plane to a radius of 1/r1/r. This does not change the magnitude frequency response and does stabilize the effect of that pole Entry 19.

A generalization of the idea of a squared-error measure is defined by raising the error to the pp power where pp is a positive integer. This error is defined by

q ( x ) = ∑ i = 0 L | H ( ω ) - H d ( ω ) | p q ( x ) = ∑ i = 0 L | H ( ω ) - H d ( ω ) | p (21)

Deczky Entry 8 developed this approach and used the Fletcher-Powell method to minimize (Equation 21). He also applied this method to the approximation of a desired group-delay function. An important characteristic of this formulation is that the solution approaches the Chebyshev or mini-max solution as p becomes large. Initial work shows the method of iteratively reweighted least squared error (IRLS) as was applied to the FIR filter design in (Reference) can also be used for LpLp and constrained least squared error optimal design of IIR filters Entry 29.

The Chebyshev Error Criterion for IIR Filter Design

The error measure that often best meets filter design specifications is the maximum error in the frequency response that occurs over a band. The filter design problem becomes the problem of minimizing the maximum error (the min-max problem).

Among several approaches to this error minimization, one is by Deczky which minimizes the p-power error of (Equation 21) for large p. Generally, p=10p=10 or greater approximates a Chebyshev result Entry 19. Another is by Dolan and Kaiser which uses a penalty-function approach.

Linear programming can be applied to this error measure by linearizing the equations in much the same way as in (Equation 15) Entry 21. In contrast to the FIR case, this can be a practical design method because the order of practical IIR filters is generally much lower than for FIR filters. A scheme called differential correction has also proven to be effective.

Although the rational approximation problem is nonlinear, an application of the Remes exchange algorithm can be implemented Entry 19. Since the zeros of the numerator of the transfer function mainly control the stopband characteristics of a filter, and the zeros of the denominator mainly control the passband, the effects of the two are somewhat uncoupled. An application of the Remes exchange algorithm, alternating between the numerator and denominator, gives an effective method for designing IIR filters with a Chebyshev error criterion. If the order of the numerator and denominator are the same and the desired filter is an ideal lowpass filter, the Remes exchange should give the same result as the elliptic function filter. However, this approach allows any order numerator or denominator to be set and any shape passband to be approximated. There are cases where a lower-order denominator than numerator results in a filter with fewer required muliplications than an elliptic-function filter Entry 19.

Prony's Method for Time-Domain Design of IIR Filters

The problem of designing an IIR digital filter with a prescribed time-domain response is addressed in this section. Most formulations of time-domain design of IIR filters result in nonlinear equations for the same reasons as for frequency-domain design. Prony, in 1790, derived a special formulation for the analysis of elastic properties of gases, which resulted in linear equations. A more general form of Prony's method can be applied to the IIR filter design by use of a matrix description Entry 2, Entry 19.

The transfer function of an IIR filter is given by

H ( z ) = B ( z ) A ( z ) = b 0 + b 1 z - 1 + ⋯ + b M z - M 1 + a 1 z - 1 + ⋯ + a N z - N = h 0 + h 1 z - 1 + h 2 z - 2 + ⋯ H ( z ) = B ( z ) A ( z ) = b 0 + b 1 z - 1 + ⋯ + b M z - M 1 + a 1 z - 1 + ⋯ + a N z - N = h 0 + h 1 z - 1 + h 2 z - 2 + ⋯ (22)

and the impulse response h(n)h(n) is related to H(z)H(z) by the zz transform.

H ( z ) = ∑ n = 0 ∞ h ( n ) z - n H ( z ) = ∑ n = 0 ∞ h ( n ) z - n (23)

Equation (Equation 23) can be written

B ( z ) = H ( z ) A ( z ) B ( z ) = H ( z ) A ( z ) (24)

which is the z-transform version of convolution. This convolution can be written as a matrix multiplication. Using the first K+1 terms of the impulse response, this is written

b 0 b 1 ⋮ b M