IEICE TRANS. FUNDAMENTALS, VOL.E89–A, NO.3 MARCH 2006
751
PAPER
Design of Equiripple Minimum Phase FIR Filters with Ripple Ratio Control∗ Masahiro OKUDA†a) , Masaaki IKEHARA† , Members, and Shin-ichi TAKAHASHI†† , Fellow
SUMMARY In this paper, we present a numerical method for the equiripple approximation of minimum phase FIR digital filters. Many methods have been proposed for the design of such filters. Many of them first design a linear phase filter whose length is twice as long, and then factorize the filter to obtain the minimum phase. Although these methods theoretically guarantee its optimality, it is difficult to control the ratio of ripples between different bands. In the conventional lowpass filter design, for example, when different weights are given for its passband and stopband, one needs to iteratively design the filter by trial and error to achieve the ratio of the weights exactly. To address this problem, we modifies wellknown Parks-McClellan algorithm and make it possible to directly control the ripple ratios. The method iteratively solves a set of linear equations with controlling the ratio of ripples. Using this method, the equiripple solutions are obtained quickly. key words: FIR filter, minimum phase, ripple ratio, bandpass filters
1.
Introduction
Since a perfect linear phase can be realized with finite impulse response (FIR) filters, they are used in many fields [1]. Efficient and optimal designs of these filters can be obtained by using an algorithm proposed by Parks and McClellan [2]. However, the linear-phase FIR filters have much redundancy, as it has linear phase also in stopbands. Thus high orders of filters are required to implement narrow transition bands, which result in large group delays. Several methods have been proposed to decrease the order and the group delay. They can be divided into two major classes. One includes the methods that try to achieve the linear phase only in the passband [3], [4]. These methods are based on approximating the complex desired response, which includes both the magnitude and the phase response. Specially, the complex Chebyshev approximation, which minimizes the Chebyshev norm of the absolute value of the complex error, is important and has been widely researched. In the second class, the phase is taken arbitrary and the focus is completely on smaller orders, thus resulting in the minimum phase FIR filters [5]–[18]. These filters are often used for Manuscript received April 4, 2005. Manuscript revised June 14, 2005. Final manuscript received November 28, 2005. † The authors are with the Faculty of Environmental Engineering, The University of Kitakyushu, Kitakyushu-shi, 808-0135 Japan. †† The author is with the Faculty of Electronics and Electrical Engineering, Keio University, Yokohama-shi, 223-8522 Japan. ∗ This paper was partially presented at IEICE Technical Report (DSP), Jan. 2004, Kansai University. a) E-mail:
[email protected] DOI: 10.1093/ietfec/e89–a.3.751
the applications, in which low delay is critical. They can be implemented by considerably smaller filter length and lower group delay as compared to the linear phase FIR filters that meet the same magnitude specification. Herrmann and Schuessler [5] propose a non-recursive design method of minimum phase FIR filters based on the Remez algorithm. They first design an optimal linear phase FIR filter for a power spectrum that is a squared ideal piecewise constant magnitude response. Then, the obtained transfer function is factored by polynomial factorization to obtain the minimum phase FIR filters. Although the designed filter using this technique [5] and some other similar techniques [6]–[8] are theoretically optimal, a serious ill-condition problem often occurs during the polynomial factorization. Several techniques [9]–[12] have been proposed to overcome this problem. Ueda et al. [9] proposed a numerical technique to robustly factor the polynomials, while Kale et al. [10] tried to facilitate the problem via balanced model truncation. Dolecek et al. [11] proposed a design based on interpolated FIR filters to reduce the problem of finding roots of higher order polynomials to lower orders. Stathaki et al. [12] presented a technique of applying Cauchy residue theorem to the logarithmic derivative of the transfer function to facilitate the factorization problem. Some alternate design techniques that avoid polynomial factorization have also been proposed for design of minimum phase filters. Boite et al. [13] and Mian et al. [14] used deconvolution of complex cepstrums. Pei et al. [15] suggested an improved design using fast Hartley transform. Venkata et al. [16] proposed a new design technique based on Hilbert transform. This method [16] is very effective in terms of complexity and accuracy. However these methods cannot directly control weight functions. For example, when one wants to obtain a filter whose passband ripples are exactly twice as large as stopband ripples, one needs to design filters repeatedly with changing weight functions. Although [17] can obtain such filters, it contains very time-consuming process. In this paper, we present a method for the design of minimum phase FIR digital filters based on controlling ripple ratios in different bands. The method iteratively solves a set of linear equations with controlling the ratio of ripples. Using this method, the equiripple solutions are obtained quickly. Our method is easily modified to the design of bandpass filters with different weights in stopband. In the next section, we briefly review a general process of the conventional minimum phase filter design methods.
c 2006 The Institute of Electronics, Information and Communication Engineers Copyright
IEICE TRANS. FUNDAMENTALS, VOL.E89–A, NO.3 MARCH 2006
752
Then, we present a new desgin method in Sects. 3 and 4. The factorization scheme is intruduced in Sect. 5. Finally, we show some examples in Sect. 6 to illustrate the advantages of our proposed method. 2.
Conventional Method
The transfer function of N-tap FIR filters is H(z) = a0 + a1 z−1 + a2 z−2 + · · · + aN−1 z−N+1 .
(1)
The filter is called a minimum phase filter if all zeros are on or inside the unit circle in Z-domain. There are many methods proposed for the design of the equiripple minimum phase filters. The Remez based algorithm is one of the most efficient method in terms of optimality and complexity. In actual, the Remez algorithm is numerically implemented by Parks-McClellan algorithm. The approach for the Remez based algorithm has three steps:
(a)
1. Design an equiripple linear phase FIR filter G(z) which has the same cut-off frequencies as H(z) and has the order of 2(N − 1) (see Fig. 1(a)). G(z) can be expressed by G(z) = b0 + b1 z−1 + · · · + b2(N−1) z−2(N−1) .
(2)
In the case that the filter is symmetry, the zero phase part of G(e jω ) is then Gm (ω) =
(b) Fig. 1
b0 + b1 cos(ω) + b2 cos(2ω)+ · · · + bN−1 cos((N − 1)ω),
(a) Gm (ω), (b) Gˆ m (ω).
where b0 = bN−1 and bn = 2bN−1−n , (n 0). ˆ by 2. As depicted in Fig. 1(b), G(z) is converted to G(z) scaling and lifting up Gm (ω) so that all the local minima of Gˆ m (ω) equal to 0 (where Gˆ m (ω) is the zero phase part ˆ jω )). of G(e ˆ to 3. Factorize G(z) ˆ = Gˆ in (z)Gˆ on (z)Gˆ out (z), G(z) where Gˆ in (z), Gˆ on (z), and Gˆ out (z) have all zeros inside, on, and outside the unit circle, respectively. Since all ˆ are 0, Gˆ on (z) has double roots on the minima of G(z) the unit circle (Fig. 2), (z − e jφi )2 . Gˆ on (z) =
Fig. 2 Conventional method: black dots, white circle denote double and single zeros.
i
Thus by picking up one of the zero pairs and the zeros of Gˆ in (z), the optimal minimum phase filter H(z) is derived as, (z − e jφi ) · Gˆ in (z) H(z) = i
3.
Design Algorithm
3.1 Design for Minimum Phase Filters with Control of Weight Functions Consider the design problem for an equiripple lowpass filter whose stopband and passband ripples are δ s and δ p , respectively. Then the ripples in stopband and passband of the ˆ jω ) are (see [9]) corresponding G(e
OKUDA et al.: MINIMUM PHASE FIR FILTERS
753
∆s =
4δ p δ2s , ∆ p = 2 ∆s, 2 2 2 + 2δ p − δ s δs
(3)
method updates the weights of G(z) according to (7). Equation (9) in the Parks-McClellan algorithm is modified to N−1
respectively. In a case that the weights of the passband and stopband are constant, we let the ratio of these be ws α= . (4) wp
n=0
δp δs
(5)
holds. Substituting (5) to (3), we get ∆s =
δ2s 4α , ∆p = ∆s. δs 2 + (2α2 − 1)δ2s
(6)
α =
∆ p 4α = . ∆s δs
(7)
Correspondingly, 2∆ s , and 1 + ∆ s (1 − 2α2 ) 2∆ s δp = α 1 + ∆ s (1 − 2α2 )
δs =
(8)
hold. When a passband and stopband ripples, δ p and δ s , are given as a specification, of course one cannot obtain the exact values of ripples. This fact and (7) imply that when designing G(z), one cannot control the ratio α, since α depends on δ s . To obtain the exact values of weights of H(e jω ), one has to design the filter iteratively with changing the weights of H(e jω ). In this paper, to address the problem, we introduce the modified Parks-McClellan algorithm. Our aim is to design the minimum phase filters with the specified α. To design ˆ G(z), the conventional Parks-McClellan algorithm solves the following set of equations iteratively with updating frequency sampling points to interpolate. N−1 n
bn cos(nωl ) +
ws (−1)l ∆ = 1, wp
bn cos(nωl ) + (−1)l ∆ = 0,
(10)
is an stopband ripples derived at the previous where δk−1 s iteration and is calculated by (8). We iteratively solve (10) until an equiripple solution is obtained. Obviously once the algorithm converges to it, the obtained filter meets the desired ratio of ripples. Bandpass Filters with Two Different Stopband Ripples
The conventional method described in Sect. 3 cannot directly be applied to the design of bandpass filters with two ˆ different stopband ripples, since G(z) obtained by lifting does not have double zeros in the stopband with smaller ripples. In this section, we modify the method in Sect. 3 to design such filters. Let the passband and the two stopband ripples of H(z) be δ p , δ s1 , δ s2 , respectively. We define δ s1 and δ s2 as the smaller and larger ripples, that is, δ s2 = βδ s1 , (β > 1). The corresponding ripples of G(z) is given by δ2s1
∆ s1 =
, 2 + 2δ2p − δ2s1 4δ p ∆ p = 2 ∆ s1 , ∆ s2 = β2 ∆ s1 . δ s1
(11)
We further define the ratio, α=
δp , δ s1
(12)
then ∆ p is expressed by ∆p =
4α ∆ s1 . δ s1
d = (β2 − 1)∆ s1
n=0
ωl ∈ stopband
ωl ∈ stopband.
(13)
Assume that we have obtained G(z) whose minima are all identical. If the desired magnitude of the stopband with smaller ripples is 0 (see Fig. 3), the one with larger ripples is
ωl ∈ passband N−1
bn cos(nωl ) + (−1)l ∆ = 0,
n=0
4.
Then the ratio of ∆ s and ∆ p is
4α ∆ = 1, δk−1 s
ωl ∈ passband N−1
where, w p and w s are weights in the passband and the stopband, respectively. When the filter is optimal, that is, perfectly equiripple, then α=
bn cos(nωl ) + (−1)l ·
(9)
in which the weights of the passband and the stopband are 1 and w p /w s , respectively. In the design of G(z), the desired weight depends on obtained ripples. Thus to control the weight of H(z), our
(14)
According to (3), the desired magnitude response depends on obtained ripples. Thus to control the local minima of H(z), our method updates the desired response of the larger ripples. The linear equations of Parks-McClellan algorithm is modified to
IEICE TRANS. FUNDAMENTALS, VOL.E89–A, NO.3 MARCH 2006
754
ˆ jω )|+ θ (ω) = − j · DFT sgn(ω) · IDFT log |G(e (17) Then the filter coefficients are derived by IDFT of ˆ jω )| + · e jθ (ω) . H(e jω ) = log |G(e
Fig. 3 N−1
This method stably solves the factorization problem [16]. Of course, one can reduce computational complexity by FFT. Finally we summarize our design algorithm for the lowpass filters: [Design Algorithm (lowpass filters)]
Bandpass filter with two different stopband ripples.
bn cos(nωl ) + (−1)l
n=0
1. Passband and stopband edges, filter length N and weights of the bands are given. And set the initial values, k = 1, and δ0s = , (: an arbitrary small number). 2. To obtain G(z), the filter of order 2(N − 1), solve the minimax problem using (10). 3. If the maximum errors are unchanged, go to next step. Otherwise, set k=k+1 and then go back to previous step. 4. Procedure 2 in the algorithm of Sect. 2 is employed. 5. The phase response is calculated by (17). 6. The filter coefficients are obtained by using IDFT of (18)
δk−1 s1 ∆ = 1, 4α
ωl ∈ passband N−1
bn cos(nωl ) + (−1)l ∆ = 0,
n=0
ωl ∈ stopband with smaller ripples. N−1
bn cos(nωl ) + (−1)l β2 ∆ = (β2 − 1)∆k−1 ,
n=0
ωl ∈ stopband with larger ripples.
(15)
where ∆k−1 is the obtained error ∆ at the previous iteration. Note that the algorithm is easily applied to any multiband filters with more than three bands. 5.
(18)
Factorization
To obtain H(z) from G(z), we use a factorization method based on the Hilbert transform, which is essentially same ˆ is strictly minimum phase (that as [16]. If a FIR filter G(z) is, zeros on the unit circle are not allowed), the magnitude ˆ jω )| and the phase response θ(ω) has a one-toresponse |G(e one corresponding given by the Hilbert transform ˆ jω )| , θ(ω) = − j · DFT sgn(ω) · IDFT log |G(e (16)
6.
Experimental Results
In this section, we show some design examples to verify the validity of our method. All examples are simulated on a Pentium IV 1.70 GHz PC. Example.1: We design a lowpass FIR filter with the specifications: 1 0 ≤ ω ≤ 0.37π , |H(e jω )| = 0 0.6π ≤ ω ≤ π N = 20 and α = 2. The values δ p and δ s of the obtained filter are 0.00253 and 0.00127, respectively. The obtained δ p /δ s is 1.99.The magnitude response and the zero plot are shown in Figs. 4 and 5. This example takes about 0.8 seconds to converge.
where
0, ω = 0, π 1, 0