Nội dung

Karam, L.J.; McClellan, J.H.; Selesnick, E.W. & Burrus, C.S. “Digital Filtering” Digital Signal Processing Handbook Ed. Vijay K. Madisetti and Douglas B. Williams Boca Raton: CRC Press LLC, 1999 c 1999 by CRC Press LLC 11 Digital Filtering 11.1 Introduction 11.2 Steps in Filter Design Lina J. Karam Arizona State University James H. McClellan Georgia Intitute of Technology Creating the Design Specifications • Specs Derived from Analog Filtering • Specifying an Error Measure • Selecting the Filter Type and Order • Designing the Filter • Realizing the Designed Filter 11.3 Classical Filter Design Methods FIR Design Methods • IIR Design Methods 11.4 Other Developments in Digital Filter Design Ivan W. Selesnick Polytechnic University C. Sidney Burrus Rice University 11.1 FIR Filter Design • IIR Filter Design 11.5 Software Tools Filter Design: Graphical User Interface (GUI) • Filter Implementation References Introduction Digital filters are widely used in processing digital signals of many diverse applications, including speech processing and data communications, image and video processing, sonar, radar, seismic and oil exploration, and consumer electronics. One class of digital filters, the linear shift-invariant (LSI) type, are the most frequently used because they are simple to analyze, design, and implement. This chapter treats the LSI case only; other filter types, such as adaptive filters, require quite different design methodologies. An LSI digital filter can be uniquely identified in the time/space domain by its impulse response h(n) (where n is an integer index). Alternatively, the LSI digital filter can be uniquely characterized in the frequency domain by its frequency response H (ω) (where ω is a real-valued frequency variable in radians), which is also the Discrete-Time Fourier Transform (DTFT) of the sequence h(n). LSI digital filters are of two main types: Finite-duration Impulse Response (FIR) filters for which the impulse response h(n) is non-zero for only a finite number of samples, and Infinite-duration Impulse Response (IIR) filters for which h(n) has an infinite number of non-zero samples. In the FIR case, the samples of the sequence h(n) are commonly referred to as the filter coefficients; for the IIR case, the filter coefficients include feedback terms in a difference equation. Digital filter design has been extensively addressed within the last 25 years. The design and realization of digital filters involve a blend of theory, applications, and technologies. For most applications, it is desirable to design frequency-selective filters which alter or pass unchanged different frequency components. In this case, the desired design specifications are given in the frequency domain by specifying a desired frequency response D(f ). Note that D(f ) is, in general, complex valued, consisting of a desired magnitude response |D(f )| and a desired phase response 6 D(f ). c 1999 by CRC Press LLC One of the most important problems is the design of a highly frequency-selective filter with sharp cutoff edges (short transition bands). However, ideal sharp edges correspond mathematically to discontinuities and cannot be realized in practice. Therefore, the filter design problem consists in finding an implementable filter whose order is low and whose frequency response H (f ) best approximates the specified ideal magnitude and phase responses which are given as the desired design specifications or constraints. The design of digital filters is typically done by performing the following steps: 1. Convert the desired design constraints into precise specifications of the desired magnitude and phase responses, designed filter type (FIR or IIR), filter order, error tolerance, or criteria. 2. Approximate the design specifications (of Step 1) by finding the implementable FIR or IIR filter such that the obtained filter frequency response best meets the design specs according to a mathematical error criterion. 3. Realize the filter using the digital technology most suitable for the considered application. While Step 2 is performed using mathematical optimization and approximation methods, Step 1 is highly dependent on the application and the detail provided by the user. Step 3 depends on the technology or software used to build the filter. Nowadays, the optimization needed in Step 2 is usually done with computer software that implements sophisticated numerical optimization routines. In addition, these design packages usually have a convenient graphical user interface to aid in the conversion of specs needed in Step 1. With such software, a filter design can be carried out quickly so that many designs can be tried in the process of getting the best filter. Since most filter design techniques involve the trade-off among competing parameters, the software can also incorporate design rules that allow the user to predict the order needed for certain specs without actually designing the filter, for example. This chapter is organized as follows. Section 11.2 provides a discussion of Steps 1 and 3, including creating the design specifications, selecting the filter type and order, specifying the error tolerances and criteria, and realizing the designed filter. Step 2 is treated in Sections 11.3 and 11.4. Section 11.3 describes the classical FIR and IIR design methods. Section 11.4 presents nonclassic and more recently developed design methods with added efficiency and/or flexibility. Finally, Section 11.5 gives examples of some of the currently available software design tools and describes the characteristics that a user can expect from such tools. 11.2 Steps in Filter Design Lina J. Karam The general filter design problem can be briefly stated as follows. Given some ideal frequency response, D(ω), find a realizable IIR or FIR digital filter whose frequency response, H (ω), approximates D(ω). The realizable filter is found by optimizing some measure of the filter’s performance, e.g., minimizing the filter order (IIR) or the filter length (FIR), or minimizing the width of the transition bands, or reducing the passband error and/or stopband error. Setting up the specifications for the general filter design problem will define these parameters and show which trade-offs are possible. 11.2.1 Creating the Design Specifications Since the frequency response of a digital filter is always periodic in the frequency variable ω with a period of 2π, the design specifications need only be specified for one period; usually, over the frequency region [−π, π]. Furthermore, when the frequency response is conjugate-symmetric (i.e., c 1999 by CRC Press LLC D ∗ (ω) = D(−ω)), then it is sufficient to specify the response only on the positive frequency interval [0, π]. The conjugate-symmetric case is the most common, because it corresponds to filters with real coefficients. The simplest case is that of an ideal low-pass digital filter with zero phase, whose frequency response can be expressed as:  1, |ω| < ωc (11.1) D(ω) = 0, ωc < |ω| < π where ωc is the cutoff frequency corresponding to the location of a sharp cutoff edge, as shown in Fig. 11.1(a). In this case, the frequency response, D(ω), is real-valued and, therefore, corresponds also to the magnitude response of the filter (since the phase is zero). Ideal frequency responses of other commonly used frequency-selective filters are shown in Fig. 11.1. FIGURE 11.1: Common ideal digital filter types. These ideal filters have frequency responses with sharp cutoff edges (discontinuities) and cannot be implemented directly. They must be approximated with a realizable system—the sharp cutoff edges need to be replaced with transition bands in which the designed frequency response would change smoothly in going from one band to the other. So, design templates need to be provided where the sharp cutoff edges are replaced with non-zero width transition bands located around the ideal cutoff edges. A typical design template for a lowpass filter is shown in Fig. 11.2, where: • ωp is the passband cutoff frequency. • ωs is the stopband cutoff frequency. The cutoff frequency ωc is usually taken to be midway between the passband and stopband cutoff frequencies. • The open interval (ωp , ωs ) is the transition band of width 1ωt = ωs − ωp . In the common design methods, no design specifications are given in the transition bands which are therefore commonly known as “don’t care bands.” However, it is usually desirable to have the frequency response change smoothly (i.e., no fluctuations or overshoots) in the c 1999 by CRC Press LLC transition bands; this requirement might not be satisfied by a design method that places no design constraints on the frequency response in the transition bands. • δp is known as the passband ripple and is the maximum allowable error in the passband. • δs is known as the stopband ripple and is the maximum allowable error in the stopband. FIGURE 11.2: Design template for a lowpass filter. The objective of filter design then is to find a realizable FIR or IIR filter whose frequency response H (ω) approximates the specified design constraints given by the design template. Ideally, the filter design process would make each of the following parameters as small as possible: δp , δs , 1ωt , IIR filter order (number of poles of H (z) which is a rational function) or FIR filter length (number of zeros of H (z) which is a finite polynomial). Practically, the filter design process minimizes one of these parameters while holding the others fixed. Traditionally, many of the filters designed in practice are specified in terms of constraints on the magnitude response and no constraints on the phase response other than those imposed implicitly by stability and/or causality requirements (e.g., poles inside unit circle in the complex Z-plane for IIR, and linear-phase for FIR [1]). More recently, design methods that include phase design specifications have been presented [2, 3, 4, 5]. In this latter case, two design templates must be provided, one for the magnitude response and another for the (passband) phase response. An ideal phase response is most likely a constant slope phase function: 6 D(ω) = −Mω The parameter M is equivalent to the desired delay of the filter (in samples). An error template for the phase would be a tolerance about the desired phase, e.g., δφ would denote the maximum allowable phase ripple, so that we require |6 H (ω) − 6 D(ω)| < δφ 11.2.2 Specs Derived from Analog Filtering Often, the desired design specifications are not given directly in the digital domain. Instead, an equivalent analog filtering operation is desired but is to be performed using an embedded digital filter. Figure 11.3 shows a standard system for processing continuous-time (-space) signals using a digital filter. The analog input signal is first transformed into a digital signal through an analog-todigital (A/D) conversion operation; then, filtering is carried out using a digital filter; finally, the filtered digital output is converted back to the analog domain using a digital-to-analog (D/A) converter. For c 1999 by CRC Press LLC this system, if the sampling period Ts of the A/D and D/A converters is chosen appropriately to avoid aliasing of the input spectrum, the overall system (consisting of the A/D converter, the digital filter, and the D/A converter) behaves as an equivalent analog filter. In this case, the frequency response Ha () of the equivalent analog filter is related to the frequency response H (ω) of the digital filter through a simple linear scaling relation between the digital frequency ω and the analog frequency . This linear scaling relation is given by (11.2) ω = Ts leading to the following expression of the analog Ha () in terms of the digital H (ω):  π    H (Ts ), || < Ts Ha () = π    0, || ≥ Ts (11.3) FIGURE 11.3: Standard system for processing analog signals using a digital (discrete-time) filter. Equivalently, H (ω) can also be expressed in terms of Ha () as follows: H (ω) = Ha (ω/Ts ), |ω| < π. (11.4) A typical filter design problem corresponding to this system is to design the digital filter such that the overall equivalent analog filter best approximates some ideal analog specifications. So, if we are given the desired analog specifications of the overall analog system, these can be turned into specifications for the desired digital filter by using Eq. (11.4). Then, a digital filter H (ω) can be designed to approximate the derived desired digital specifications. Finally, the resulting analog frequency response of the overall system can be found using Eq. (11.3), for example, to compare with the ideal analog response. 11.2.3 Specifying an Error Measure An error measure is needed to assess how much the designed filter H (ω) deviates from the desired filter D(ω). Defining the pointwise error E(ω) as: E(ω) = [D(ω) − H (ω)], (11.5) we must reduce E(ω) to a scalar error measure (also called an error norm). With a correctly chosen norm, there are many possible optimization algorithms that will compute the best filter parameters to minimize the chosen error norm. The following error norms are the most commonly used in filter design: c 1999 by CRC Press LLC • Mean Squared Error (MSE) or L2 norm  1/2 Z 1 2 |E(ω)| dω E2 = 2π B (11.6) • Lp norm which is a generalization of the L2 norm and where p is a non-zero integer  1 Ep = 2π • Chebyshev or L∞ norm Z p B 1/p |E(ω)| dω E∞ = max |E(ω)| ω∈B (11.7) (11.8) The Chebyshev error norm limits the worst case deviation from the ideal specifications. In the above definitions, | · | denotes the complex error magnitude and B is the frequency region of interest over which the error norm is to be minimized. The frequency subset B ⊂ [−π, π ) is taken to be the union of the desired passbands and stopbands. A more selective control of the approximation accuracy can be achieved by introducing a weighting function W (ω) in Eq. (11.5) as follows: E(ω) = W (ω)[D(ω) − H (ω)]. (11.9) The weighting function W (ω) must be a real, strictly positive and continuous function on B. It can force a better match over selected regions or frequency points relative to other regions in B. Alternatively, note that Eq. (11.5) reduces to Eq. (11.9) if we replace D(ω) with W (ω)D(ω) and H (ω) with W (ω)H (ω). 11.2.4 Selecting the Filter Type and Order As mentioned in Section 11.1, there are two main types of filters, namely FIR and IIR. These differ in their characteristics and in the way they are designed. Since the design algorithm depends strongly on the choice of IIR vs. FIR filter, the designer should make this decision as early as possible. Although the desired frequency response specifications can be approximated with either type of filter, deciding which of the two filter types to use depends on many factors including the implementation hardware, as well as the magnitude and phase characteristics of the resulting filter. To aid in this decision, the main characteristics of FIR and IIR filters are discussed below. 11.2.4.1 FIR Characteristics 1. The impulse response h(n) has a finite length, i.e., h(n) is non-zero only for a finite range of indices n. For a general N-length FIR system, h(n) 6 = 0 only for N1 ≤ n ≤ N2 = (N1 + N − 1). When N1 ≥ 0, the filter is also causal. 2. The FIR frequency response H (ω) is a finite-degree polynomial in ej ω of the form H (ω) = N2 X h(n)(ej ω )−n (11.10) n=N1 where N1 and N2 are (negative or positive) integers corresponding to the indices of the first and last samples of h(n), respectively. The N impulse response samples are the free parameters of the design procedure. This form is general enough to represent non-causal filters such as zero-phase filters. c 1999 by CRC Press LLC 3. Designing an FIR filter consists in finding the polynomial H (ω) that best approximates the design specifications. This is done by computing the “optimal” (relative to some criteria) 2 impulse response samples {h(n)}N n=N1 , which correspond to the unknown coefficients of the polynomial H (ω). The impulse response length N is usually fixed, but it could also be considered as a free parameter to be optimized. Procedures for designing FIR filters are given in Sections 11.3.1 and 11.4.1. 4. The filter transfer function, denoted by H (z), is the z-transform of h(n) and is useful for studying the stability of the system. For FIR filters, H (z) is a finite-degree polynomial in the complex variable z and is given by jω H (z) = H (e )|ej ω =z = 5. 6. 7. 8. 9. N2 X h(n)z−n . (11.11) n=N1 It follows that the function H (z) has no poles except possibly at 0 or ∞, i.e., it cannot be infinite for any point z with 0 < |z| < ∞. It has only zeros (points z at which H (z) = 0). Therefore, an FIR filter is always stable. FIR filters allow the design of causal linear-phase systems which are very important and widely used in practice. In fact, in many signal processing applications, such as speech and image processing, it is desirable to pass some portion of the signal frequency band with minimal distortion. For that purpose, linear-phase systems are particularly desirable since the effect of the linear-phase is a pure time delay. For a more detailed discussion of linear-phase systems, the reader is referred to [1]. Because the impulse response is of finite length, FIR filters are realized using the convolution operation [1] which can be implemented directly in the time/space domain, or in terms of the FFT in the frequency domain. More details about the implementation will be given in Section 11.2.6. Since FIR filters have no feedback loops, they are relatively insensitive to round-off noise. Noise due to coefficient quantization can be a problem for very long filters, but can be mitigated by avoiding the direct-form structures, and using special structures such as the cascade form for implementation. FIR filters with very long impulse responses (N ≈ 500) might be required to meet certain design specifications, e.g., high accuracy and/or short transition bands. Longer filters lead to an increased complexity for both design and implementation. They require significant computing time to optimize all the parameters h(n), and also many operations per second in the actual filter implementation. The trade-off among the filter design parameters has been determined empirically for some types of FIR designs. The following simple (approximate) formula shows the relationship among the ripples, bandedges, and filter length (N ) for one method, the Parks-McClellan algorithm : p −20 log10 δp δs − 13 (N − 1)1ω ≈ 2.324 where 1ω = ωs − ωp is the transition width. This formula allows the designer to predict the value of N that will be needed to satisfy specs given for { ωp , ωs , δp , δs }. Other design formulas are given in Section 11.3.1. c 1999 by CRC Press LLC 11.2.4.2 IIR Characteristics 1. The impulse response h(n) has an infinite number of non-zero samples (infinite length). As an example, for a general IIR filter, h(n) 6 = 0 only for No ≤ n ≤ ∞, where No is a non-negative integer (commonly, No is taken to be 0; in this case, the filter is said to be causal). 2. The frequency response H (ω) is a rational function, i.e., a ratio of two finite-degree polynomials in ej ω of the form H (ω) = PM bk e−j ωk B(ω) = e−j ωNo Pk=0 N −j ωk A(ω) k=0 ak e (11.12) where No is an integer constant. The order of an IIR filter is equal to N, which is the degree of the denominator in Eq. (11.12); usually the degree of the numerator M is no greater than N. The order N also determines the number of previous output samples that need to be stored and then fed back to compute the current output sample. Therefore, IIR systems are also known as feedback systems. The filter coefficients {bn } and {an } in Eq. (11.12) correspond to the unknown (free) parameters of the design. 3. Designing an IIR filter amounts to finding the rational function H (ω) that best approximates the design specifications. In the frequency domain, this is done by computing the “optimal” (relative to some criteria) coefficients {bn } and {an } in Eq. (11.12) for the rational function H (ω). The filter order N is usually fixed, but can also be considered as a free parameter to be optimized. Procedures for designing IIR filters are given in Sections 11.3.2 and 11.4.2. 4. As mentioned previously, the filter transfer function, denoted by H (z), is the z-transform of h(n) and is useful for studying the stability of the system. In the context of LSI filters, stability implies that a bounded input to the filter will always result in a bounded output. For IIR filters, H (z) is a rational function in the complex variable z and is given by PM   bk z−k H (z) = H ej ω |ej ω =z = z−No Pk=0 N −k k=0 ak z (11.13) The roots of the denominator polynomial are poles of the function H (z), i.e., H (z) is infinite at points z with 0 ≤ |z| < ∞. Stability then requires that no poles lie on the Unit Circle (U.C.) (|z| = 1) in the z-plane. Causality and stability require that the poles lie inside the U.C. in the z-plane. So, it is possible to obtain a resulting IIR filter that is unstable. Also, coefficient quantization noise might severely affect the response of the filter and its stability by disturbing the poles locations and by driving some of the poles closer to or onto the U.C. 5. It is not possible to design causal linear-phase IIR filters. The resulting IIR causal realizable filters must have a non-linear phase response. Forward-backward filtering can be used as an implementation to approximate a zero-phase response [1]. 6. Because the impulse response is infinitely long, convolution can no longer be used to implement the IIR filters. Instead, IIR filters are efficiently implemented using feedback difference equations as described in Section 11.2.6. 7. The noise characteristics of an IIR filter can be a major consideration when doing an implementation, especially in fixed-point arithmetic. Coefficient quantization degrades the actual filter response from that designed by high-precision software. More critical is c 1999 by CRC Press LLC round-off noise sensitivity which can be amplified by the feedback loops in the filter. 8. Compared to FIR filters, IIR filters can achieve the desired design specifications with a relatively low order (as few as 4 to 6 poles). So, fewer unknown parameters need to be computed and stored, which might lead to a lower design and implementation complexity. However, the phase response of IIR filters is never linear, which leads to the use of allpass filters to compensate the group delay, and thus raises the order of the filter and the complexity of the design process. 9. IIR filters are commonly designed by using closed-form design formulas corresponding to classical filter types. While for FIR filters the length-estimating formulas are only approximate, the order-estimating formulas for IIR filters are exact since they are derived from the mathematical properties of the classical prototypes. These formulas are very useful to obtain the IIR filter order needed to satisfy the desired design specifications. 11.2.5 Designing the Filter After the designed filter type (FIR or IIR) is specified, a suitable design procedure can be selected depending on the chosen filter type. Popular design procedures are based on computing the unknown filter parameters by optimizing one of the error criteria indicated in Section 11.2.3. For FIR filters, the two main classical methods are the windowing method [1] and the ParksMcClellan (Remez) algorithm [6]. The windowing method minimizes the MSE when a rectangular window (corresponding to pure truncation of the ideal impulse response) is used at the expense of possible large overshoots near the band edges and large ripples in the resulting frequency response. It is suboptimal when other general windows are used. However, the edge overshoot, transition width, and ripple height can be controlled by using different types of windows as described in Section on page 11-13. The Parks-McClellan (Remez) algorithm minimizes the Chebyshev (L∞ ) error norm resulting in optimal equiripple designs. However, the original Parks-McClellan algorithm is restricted to the design of linear-phase filters with a symmetric magnitude response. An extension of this algorithm that allows the design of optimal FIR filters with arbitrary magnitude and phase specifications has been presented by Karam and McClellan in [2, 3]. Linear-programming-based [4, 7] and Constrained least square [8] optimization methods also have been presented to allow the inclusion of additional important design constraints. These and other FIR design procedures are described in Sections 11.3.1 and 11.4.1. While the design of FIR filters is typically performed directly in the digital domain, IIR filters are commonly designed by transforming the digital design specifications into analog design specifications and performing the filter design in the analog domain. The resulting analog filter is then transformed into a digital filter using a suitable transformation. One important classical IIR design method is the Bilinear Transformation method. Digital-only IIR design methods have also been presented. A description of IIR design procedures is given in Sections 11.3.2 and 11.4.2. 11.2.6 Realizing the Designed Filter Realizing the designed digital filter corresponds to computing the output of the filter in response to any given input. For LSI filters, this is simplified by the fact that the input and output signals are related through a simple convolution operation in the time/space domain. If x(n) is the input, y(n) the corresponding output, and h(n) the impulse response of the LSI filter, then this relation is given by y(n) = h(n) ∗ x(n) = N2 X k=N1 c 1999 by CRC Press LLC h(k)x(n − k), (11.14)