Design of Optimum High Order Finite Wordlength Digital FIR Filters ...
EURASIP SIGNAL PROCESSING
1
Design of Optimum High Order Finite Wordlength Digital FIR Filters with Linear Phase Gennaro Evangelista Abstract— A novel iterative quantization procedure for the design of finite wordlength linear phase FIR filters of high order and minimum frequency domain error is proposed: a one by one increased number of filter coefficients is quantized where the augmented frequency domain error is re-minimized in each case. This new approach achieves a smaller frequency domain error than rounding of the optimal non-quantized coefficients. It is also successfully applied to the design of Lth band filters with minimum frequency domain error. FIR filters up to a filter order 1500 are designed. Keywords— Design of high order FIR filters, Constrained Least Squares (CLS) design, finite wordlength coefficients
I. Introduction Conventional design methods of FIR filters with linear phase minimize a norm of the frequency domain error E(Ω) = H(Ω) − D(Ω)
(1)
on a prescribed approximation domain B, being a subset of [0, π]. H(Ω) is the real-valued amplitude response and D(Ω) a given real-valued target function. (In this contribution no weighting function is considered.) Two error norms are commonly used: the L∞ -error norm (maximum error) ||E(Ω)||∞ = max |E(Ω)| Ω∈B
and the L2 -error norm (mean squared error) Z 1 2 ||E(Ω)||2 = |E(Ω)|2 dΩ. 2π Ω∈B
(2)
(3)
The coefficients of a filter with minimum L∞ -error norm are computed by the well known MPR-program [16]. The filter coefficients for a minimum L2 -error norm are obtained by solving a system of linear equations [3]. Both algorithms yield infinite wordlength coefficients. However the realization of a digital filter with finite expenditure allows inevitably only filter coefficients of finite precision. If a digital signal processor with floating point arithmetic is used, the effects of a finite coefficient wordlength w can (nearly always) be neglected. However for a fast and cheap implementation of a digital filter (e.g. as ASIC) a fixed point arithmetic has to be used (as assumed in the following) and the effects of finite coefficient wordlength must be considered. This is especially true for very long filters with applications to sharp cut-off or sample rate alteration [24]. The author is with SIEMENS ICM, Munich. email: [email protected]. This work was performed during his time as research assistant with Digital Signal Processing Group at the Ruhr-Universit¨ at Bochum, Germany and was supported by Deutsche Forschungsgemeinschaft under contract GO 849/1-1.
The most widespread method of finding finite wordlength coefficients is the direct quantization (DQ) method in which the infinite wordlength coefficients, obtained by one of the aforementioned filter design methods [16], [19], are rounded to yield quantized coefficients. As a result the respective error norm increases, possibly beyond the maximally allowed value. To achieve a smaller error norm the coefficient wordlength and/or the filter order must be increased [9], leading to an additional expenditure. As the finite wordlength coefficients found by the direct quantization method are not optimal in sense of the error norm (since this method does not take into account the coefficient wordlength), finite wordlength filters of original order and coefficient wordlength may exist with smaller error norm. Several contributions deal with the problem of finding the optimal finite wordlength coefficients of FIR filters with minimum L∞ -error norm by using mixed integer linear programming [10], [14], [15], simulated annealing [2] and optimal or local search methods [4], [11], [25]. All these methods have nothing in common with the conventional infinite wordlength FIR filter design method [16] and are very time-consuming or, even worse, do not converge for high filter orders. A more promising approach, where the coefficients are quantized one by one and the L2 -error norm is reminimized each time, is presented in [13]. An application of this approach to the design of FIR filters with minimum L∞ -error norm requires an algorithm for re-minimization of the L∞ -error norm. A modification of the MPR algorithm [16] can not be recommended because the MPR algorithm can only influence as many extrema as free variables (number of non-quantized filter coefficients) exist, and the L∞ error norm increases significantly with a decreasing number of extrema [6], [5]. This contribution focuses on the design of linear-phase FIR filters of high order and minimum L∞ -error norm (2). A filter with minimum L∞ -error norm is further on referred to as optimal. Instead of a modified MPR algorithm a modified Constrained Least Square (CLS)-algorithm [1] is proposed for re-minimization. The filter coefficients are assumed to be identical with the impulse response h(k) as it is the case for a direct form implementation. The coefficients are restricted to take only discrete values on the quantization grid [−1 + q, −1 + 2q, . . . , 1 − q, 1] with q = 2−w+1 . First a procedure similar to the method in [13] is outlined, where the filter coefficients are iteratively quantized. Then the modified CLS-algorithm is presented which is used for the re-minimization of the L∞ -error norm after the quantization of one somehow selected coefficient. The
2
EURASIP SIGNAL PROCESSING
design examples show the power of this new approach.
error: min |∆h(k)| = min |hQ (k) − h(k)|.
II. Iterative Quantization Procedure The proposed Iterative Quantization (IQ) consists of 3 steps. Steps 2 and 3 are repeated iteratively. In each iteration one additional coefficient is quantized. So the set of quantized coefficients increases from iteration to iteration until all coefficients are quantized. The IQ procedure is outlined subsequently (Fig. 1): 1. Initialization: The L∞ -error norm is minimized with the MPR [16] yielding the non-quantized coefficient set {h0 (k)}. The set {h0Q (kf )} of quantized coefficients and the set {kf } of indices of quantized coefficients are void. 2. Selection and Quantization of an Additional Coefficient: Out of the actual set {hi−1 (k)} of non-quantized coefficients one coefficient h(kS ) is selected, rounded to the next quantization step and assigned to the next set {hiQ (kf )} of quantized (fixed) coefficients. The index kS is shifted from the set {k} of indices of non-quantized coefficients to the set {kf } of indices of quantized coefficients. Hence, from one iteration (i−1) to the next i, the set {hiQ (kf )} increases by one coefficient, while the set {hi (k)} diminishes by one coefficient. 3. Re-Minimization of L∞ -Error Norm: The L∞ -error norm increased due to the preceding coefficient quantization is re-minimized. Note that only the set {hi (k)} can be used for re-minimization. The iterative procedure ends as soon as the set {hi (k)} is void.
∀k
III. Re-Minimization of L∞ -Error Norm For the re-minimization of the L∞ -error norm considering the set {hiQ (kf )} of fixed filter coefficients a modified CLS-algorithm is used. The CLS-algorithm [1] minimizes the L2 -error norm (3) while the L∞ -error norm (2) is constrained not to exceed prescribed bounds. To use the CLSalgorithm for the minimization of the L∞ -error norm these prescribed bounds must be chosen as tighten as possible. The combination of the Iterative Quantization with the modified CLS-algorithm for the re-minimization of the L∞ error norm is subsequently called IQCLS. First the original CLS-algorithm [1] is described and then its modification. A. CLS-Algorithm for the Design of FIR Filters The amplitude response H(Ω) of a zero-phase FIR filter of even order n and even symmetry1 is H(Ω) = h(0) +
h
= [1, 2 cos(Ω), . . . , cos( n2 Ω)], T
= [h(0), . . . , h( n2 )] .
(7)
(8)
Since E(Ω) is real and a scalar:
|E(Ω)|2 = E 2 (Ω) = E(Ω)T E(Ω).
Quantization of h(kS) hQ(kS) {hQi(kf)}
i=i+1
Re-Minimization of L -error {hi(k)}
where
Fig. 1. Flow chart of the Iterative Quantization (IQ) Procedure
There exists a large number of possible selection criteria for step 2. The most successful selection strategy [6], [5] chooses the coefficient with the smallest absolute rounding
(9)
Introducing (9) with (8) in (3) leads to an alternative expression of the L2 -error norm [3] Z 1 ||E(Ω)||22 = hT Ah − 2hT b + D2 (Ω)dΩ, (10) 2π Ω∈B
(6)
(5)
With (5) the error E(Ω) of (1) is
E(Ω) = aT h − D(Ω).
2 cos(kΩ)h(k) = aT h
with2
i=0
Selection of one h(kS) out of {hi-1(k)}: kS {kf}, {k}={k} kS
n/2 X k=1
i=1
(4)
The IQ method is suboptimal but leads with a great probability to a smaller L∞ -error norm than the direct quantization of the coefficients from the MPR-program (DQMPR), because the L∞ -error norm is constantly reminimized considering the already quantized coefficients.
aT Initialization: Minimization of L - error {h0(k)} with { kf},{hQ0(kf)} void
∀k
A
=
b
=
Z 1 aaT dΩ, 2π Ω∈B Z 1 aD(Ω)dΩ. 2π Ω∈B
(11) (12)
1 Equivalent equations can easily be found for all other types of linear phase FIR-filters [19]. 2 Boldface capital letters are used for matrices, boldface lower case letters for vectors.
EVANGELISTA: DESIGN OF DIGITAL FIR FILTERS
3
The minimization of (10) with respect to h leads to the starting solution of the CLS-algorithm: hstart = A
−1
b.
(13)
If all local extrema of the amplitude response of hstart do not exceed the prescribed upper and lower bounds Rlo (Ω) and Rup (Ω) on B, no iteration is necessary and hstart is the final solution of the CLS-algorithm. Otherwise the amplitude response is set to the upper or lower bound, respectively, at all extremal frequencies Ων , where H(Ω) violates these bounds: ½ Rup (Ων ), if H(Ων ) > Rup (Ων ), ! H(Ων ) = (14) Rlo (Ων ), if H(Ων ) < Rlo (Ων ), or in matrix notation Bh = r
(15)
with
.. . T B= a (Ων ) , .. .
.. .
r= Rup,lo (Ων ) . .. .
(16)
Hence, in this case the extended L2 -error norm, the Lagrange-function [8] ELag = ||E(Ω)||22 − mT (Bh − r)
(17)
with the Lagrange-multipliers m, has to be minimized. The minimization of the Lagrange-function (17) in conjunction with (10) with respect to h leads to Ah − 0.5B T m = b.
(18)
Eqs. (15) and (18) commonly represent a set of linear equations with the solution [1]: m = 2(BA−1 B T )−1 (r − BA−1 b) h = A−1 (b + 0.5B T m).
(19) (20)
The amplitude response of h is calculated and if its extrema are beyond the bounds at the (new) extremal frequencies Ων , the procedure starting at eq. (14) is repeated. The CLS-algorithm stops if all extrema of the amplitude response are within the bounds (convergence) or the set {Ων } of extremal frequencies does not change anymore (no convergence, because the bounds are to tight). B. Modification of the CLS-Algorithm To consider fixed coefficients hQ (kf ) in the CLSalgorithm the amplitude response (5) is expressed depending of a vector hf containing the fixed coefficients hQ (kf ) and a vector hu containing all non-quantized coefficients to be used for re-minimization: T H(Ω) = aT u hu + af hf .
(21)
Replacing H(Ω) with (21) in (1) yields E(Ω) = aT u hu − DMod (Ω),
(22)
DMod (Ω) = D(Ω) − aT f hf .
(23)
where As the expressions of the modified error (22) and the original error (1) and (5) are the same replacing a with au , h with hu and D(Ω) with DMod (Ω), the modified and original CLS-algorithm are also the same by using the same replacements. So A (11), b (12) and B (16) must be replaced with Z 1 (24) au aT AMod = u dΩ, 2π Ω∈B Z 1 au DMod (Ω)dΩ, (25) bMod = 2π Ω∈B .. T. BMod = (26) au (Ων ) . .. . The matrix AMod to be inverted in the modified CLSalgorithm is of lower order than the matrix A of the original CLS-algorithm. So the consideration of the fixed coefficients hf leads to a saving in computational expenditure. The upper and lower bounds are introduced as follows Rup (Ω) = Rlo (Ω) =
D(Ω) + δ, D(Ω) − δ,
(27) (28)
with δ as input parameter of the modified CLS-algorithm. δ must be iteratively decreased for the re-minimization of the L∞ -error norm (step 3 of the IQ) until the modified CLS-algorithm does not converge anymore. The smallest value of δ with convergence of the modified CLS-algorithm is the minimal L∞ -error norm. A major advantage of the IQCLS algorithm is that is based on a conventional design method (the CLS-algorithm [1]): All experience, especially with the design of filters of high order, can be applied. IV. Examples This section focuses on the design of lowpass filters, because lowpass filters are widely used and similar results are also obtained for all other types of frequency selective filters. A first target function D(Ω) is ½ L, for 0 ≤ Ω ≤ ΩP , D(Ω) = (29) 0, for ΩS ≤ Ω ≤ π, where the factor L in D(Ω) is chosen in such a way that the quantization grid [−1 + q, −1 + 2q, . . . , 1 − q, 1] for the coefficients hQ (kf ) is optimally exploited. In Fig. 2 the magnitude resonse |H(Ω)/L| of a lowpass filter designed with IQCLS is depicted. Although already 7 of 11 coefficients are quantized with 8 bit, |H(Ω)/L| shows nearly equi-ripple behaviour.
4
EURASIP SIGNAL PROCESSING
3) outperforms the DQMPR approach, if the difference between the minimax errors of the DQMPR and MPR is large (see figures 4, 5, 6).
Magnitude / dB
0.2
0.1
0
−0.1
−10 −0.2
0
0.05
0.1
0.15
0.2
Ω/π
0.25
0.3
0.35
0.4
−15
L −error norm / dB
−35
−40
−45
−50 0.5
0.6
0.7
0.8
0.9
1
Ω/π
1.1
1.2
1.3
1.4
1.5
−20
−25
∞
Magnitude / dB
−30
−30
Fig. 2. 20 log10 |H(Ω)/L| with IQCLS for n = 20, ΩP = 0.4π, ΩS = 0.6π, 7 of 11 coefficients quantized with w = 8: passband (above), stopband (below)
−35 50
100
150
200
Filter order n
A. Design of FIR Lowpass Filters FIR lowpass filters were designed for different specifications (29) and wordlengths w using the IQCLS algorithm. The following observations can be reported: −13
n=50
• For a given wordlength w and specification D(Ω) a filter order nmax exists beyond which the L∞ -error norm does not improve anymore for both methods, IQCLS and DQPMR (see figure 5: nmax ≈ 250 for w = 8, figure 6: nmax ≈ 180 for w = 6). This result is consisting with [9].
n=60
−15
−10
−16
n=70 −20
−17
w=8 −30
n=80
−18
−40 −19
L∞−error norm / dB
∞
L −error norm / dB
−14
Fig. 4. Log. L∞ -error norm with IQCLS (◦), DQMPR (∗) and MPR (+) for ΩP = 0.192π, ΩS = 0.208π, L = 5, w = 8
n=90
−20
n=100 −21 0
5
10
15
20
25
30
35
40
45
50
Number of quantized coefficients
−60
w=16
−70 −80 −90
−100
Fig. 3. Remaining log. L∞ -error norm after IQCLS re-minimization for ΩP = 0.192π, ΩS = 0.208π, L = 5, w = 8
−110 −120
For a small number of quantized coefficients the remaining L∞ -error norm increases scarcely after re-minimization in step 3 of the IQCLS algorithm (s. figure 3). This is due to the fact that for IQCLS re-minimization enough nonquantized coefficients are left (e.g. for n = 90, 100 in figure 3). In the opposite case a smaller final L∞ -error norm can sometimes be achieved by quantizing all remaining nonquantized coefficients directly. The number of coefficients to be quantized in one step can be deduced of the L∞ -error norm of figure 3 (e.g. for n = 50, 70, 80 as indicated by the dashed lines). • The overall increase of the L∞ -error norm during the IQCLS is smaller for larger coefficient wordlength w. • As expected the IQCLS algorithm (with adequate direct quantization of the remaining coefficients as in figure
w=12
−50
•
100
200
300
400
500
600
700
800
900
Filter order n
Fig. 5. Log. L∞ -error norm with IQCLS (◦), DQMPR (∗) and MPR (+) for ΩP = 0.192π, ΩS = 0.208π, L = 5
FIR filters were successfully designed using the IQCLS algorithm (implemented in Matlab) up to a filter order of n = 1500. Beyond this upper bound for n the Choleskydecomposition [7], [8] used to invert AMod does not work properly. •
The IQCLS method needs considerably more computation time than the DQMPR (e.g. 750 re-minimizations of L∞ error are required for n = 1500 in IQCLS while in DQMPR there is only one minimization of L∞ -error).
EVANGELISTA: DESIGN OF DIGITAL FIR FILTERS
5
−10
−29.5
−30
L −error norm / dB
−30.5
−20
−25
−31
−31.5
−32
−32.5
∞
∞
L −error norm / dB
−15
−30
−33
−33.5
−34
−35 50
100
150
−34.5
200
Filter order n
The IQCLS algorithm can also be applied to the design of Lth-band (Nyquist-)filters [17] with minimum L∞ -error norm. The zero phase impulse response h(k) of Lth-band filters must hold ½ 1, for ν = 0, h(νL) = (30) 0, for ν ∈ Z{0}. Lth-band filters are advantageously used as anti-imaging filters in interpolators with upsampling factors L [24]. An appropriate target function D(Ω) for Lth-band filters is given by [22]: L, 0 ≤ Ω ≤ ΩP , 0, 2π L ν − Ωp ≤ Ω ≤
4
5
6
7
8
9
10
Fig. 7. Log. L∞ -error norm with IQCLS, (◦), CMPR (∗) and MPR π (+) for ΩP = 0.96 L , n = 40L, non-quantized coefficients
V. Conclusion and Outlook
B. Design of Lth-Band Filters
½
3
Upsampling factor L
Fig. 6. Log. L∞ -error norm with IQCLS (◦), DQMPR (∗) and MPR (+) for ΩP = 0.192π, ΩS = 0.208π, L = 5, w = 6
D(Ω) =
2
2π L ν
+ Ωp ,
(31)
with ν = 1, . . . , L−1. Obviously D(Ω) has L−1 stopbands. The task of designing Lth-band filters with a minimum L∞ -error norm is addressed by [2], [12], [15], [20], [21] using methods with high computational expenditure. The use of the IQCLS algorithm makes the design very fast: The coefficients given by (30) are all treated as fixed coefficients in step 2 and the re-minimization of the L∞ -error norm by the modified CLS algorithm is done once in step 3 yielding the optimal, non-quantized coefficients. Of course, the design procedure can be continued in order to have all coefficients quantized, as treated above. The L∞ -error norm of Lth-band filters for different upsampling factors L and design methods is depicted in figure 7. The coefficients h(k) resulting from the MPR-program with D(Ω) according to (31) violate (30). Hence the result is not an Lth band filter. The corresponding L∞ -error norm is only depicted as lower bound of the achievable L∞ error norm. The forced correction of the coefficients h(k) from the MPR to hold (30) is done by the CMPR. As can clearly be seen the IQCLS algorithm leads to a considerably smaller L∞ -error norm than the CMPR [6].
The design of linear phase FIR filters with finite wordlength coefficients and minimum L∞ -error norm using an Iterative Quantization approach combined with a modified CLS algorithm (IQCLS) has been presented. The IQCLS algorithm leads to a smaller L∞ -error norm than the DQMPR (direct quantization of the coefficients h(k) resulting from the MPR program) for the case, where a great difference between the L∞ -error norms of MPR and DQMPR exists. The IQCLS has successfully been applied to filters up to a filter order of n = 1500. It is also very suitable for the design of Lth-band filters with minimum L∞ -error norm. For a further reduction of the L∞ -error norm the direct quantization of one or more remaining coefficients at the end of the IQCLS (see sect. II) might be replaced by applying local search algorithms [11]. The IQCLS can also be easily applied to the design of filters with different coefficient wordlengths. This leads to simple and fast multiplications for those coefficients with very small wordlengths. Especially the coefficients quantized in the first iteration steps of the IQCLS require only a small wordlength [6] (see also figure 3). The design of filter with solely power-of-two coefficients or with coefficients in CSD-code representation [18] is possible, too. With a small modification, the IQCLS can also be used for the design of finite wordlength linear phase FIR filters and Lth-band filters with minimum L2 -error norm. Hence, with the IQCLS a very flexible and powerful filter design tool is given. Because the selection criterion (4) used in the second step of the IQ is not proven to be optimal, further research could deal with finding other more appropriate selection criteria. Acknowledgement The author wishes to express his appreciation to Prof. Dr.-Ing. H. G. G¨ockler and the unknown reviewers for their careful review of the manuscript and valuable suggestions.
6
EURASIP SIGNAL PROCESSING
References [1] [2]
[3] [4]
[5] [6] [7] [8] [9] [10]
[11]
[12]
[13] [14] [15] [16]
[17] [18] [19] [20] [21] [22] [23]
[24] [25]
J. W. Adams, FIR Digital Filters with Least Squares Stopbands Subject to Peak-Gain Constraints, IEEE Transactions on Circuits and Systems, 39:376-388, April 1991. N. Benvenuto, M. Marchesi, and Aurelio Uncini, Applications of Simulated Annealing for the Design of Special Digitals Filters, IEEE Transactions on Signal Processing, 40:323-332, February 1992. C. S. Burrus, A. W. Soewito, R. A. Gopinath, Least Squared Error FIR Filter Design with Transition Bands, IEEE Transactions on Signal Processing, 40:1327-1340, June 1992. C.-L. Chen and A. N. Willson Jr., A Trellis Search Algorithm for the Design of FIR Filters with Signed-Powers-of-Two Coefficients, IEEE Transactions on Circuits and Systems II, 46:2939, January 1999. M. Erdogan, Vergleichende Untersuchung zum Entwurf linearphasiger FIR-Filter mit quantisierten Koeffizienten, Diplomarbeit, Ruhr-Universit¨ at Bochum, 1999. G. Evangelista, Zum Entwurf digitaler Systeme zur asynchronen Abtastratenumsetzung, Ph.D. Thesis, Bochum, 2000. F. R. Gantmacher, Matrizentheorie, Berlin: Springer, 1986. G. H¨ ammmerlin and K.-H. Hoffman, Numerische Mathematik, Berlin: Springer, 1991. ¨ U. Heute, Uber Realisierungsprobleme bei nicht-rekursiven Digitalfiltern, Ph.D. Thesis, Erlangen, 1975. D. M. Kodek, Design of Optimal Finite Wordlength FIR Digital Filters using Integer Programming Techniques, IEEE Transactions on Acoustics, Speech, Signal Processing, 28:304-308, June 1980. D. M. Kodek and K. Steiglitz, Comparison of Optimal and Local Search Methods for Designing Finite Wordlength FIR Digital Filters, IEEE Transactions on Circuits and Systems, 28:28-32, January 1981. J. K. Liang, R. J. P. DeFigueiredo, and F. C. Lu, Design of Optimal Nyquist, Partial Response, N th Band, and Nonuniform Tap Spacing FIR Filters using Linear Programming Techniques, IEEE Transactions on Circuits and Systems, 32:386392, April 1985. Y. C. Lim and S. R. Parker, A Discrete Coefficient FIR Digital Filter Design Based Upon an LMS Criteria, Int. Symposium Circuits and Systems, Rome, Italy, 3:796-799, May 1982. Y. C. Lim and S. R. Parker, FIR Filter Design Over a Discrete Powers-of-Two Coefficient Space, IEEE Transactions on Acoustics, Speech, Signal Processing, 31:583-591, June 1983. Y. C. Lim and B. Liu, Design of Cascade FIR Filters with Discrete Valued Coefficients, IEEE Transactions on Acoustics, Speech, Signal Processing, 36:1735-1739, November 1988. J. H. McClellan, T. W. Parks and L. R. Rabiner, A Computer Program for Designing Optimum FIR Linear Phase Digital Filters, IEEE Transactions on Audio Electroacoustics, 21:506-526, 1973. F. Mintzer, On Half-band, Third-Band, and N th Band FIRFilters and Their Design, IEEE Trans. on Acoustics, Speech, Signal Processing, ASSP-30:734-738, October 1982. S. K. Mitra and J. F. Kaiser, Handbook for Digital Signal Processing, p. 206, J. Wiley & Sons: New York, 1993. T. W. Parks and C. S. Burrus, Digital Filter Design, New York: John Wiley & Sons, 1987. H. Samueli, On the Design of Optimal Equiripple FIR Digital Filters for Data Transmission Applications, IEEE Transactions on Circuits and Systems, 35:1542-1546, December 1988. T. Saram¨ aki and Y. Neuvo, A Class of FIR Nyquist (N thBand) Filters with Zero Intersymbol Interference, IEEE Transactions on Circuits and Systems, 34:1182-1190, October 1987. R. W. Schafer and L. R. Rabiner, A Digital Signal Processing Approach to Interpolation, Proceedings of the IEEE, 61:692702, June 1973. I. W. Selesnick, M. Lang and C. S. Burrus, Constrained Least Square Design of FIR Filters without Specified Transition Bands, IEEE Transactions on Signal Processing, 44:1879-1892, August 1996. P. P. Vaidyanathan, Multirate Systems and Filter Banks, Englewood Cliffs: Prentice Hall, 1993. Q. Zhao and Y. Tadakoro, A Simple Design of FIR Filters with Powers-of-Two Coefficients, IEEE Transactions on Circuits and Systems, 35:566-570, May 1988.
1
Design of Optimum High Order Finite Wordlength Digital FIR Filters with Linear Phase Gennaro Evangelista Abstract— A novel iterative quantization procedure for the design of finite wordlength linear phase FIR filters of high order and minimum frequency domain error is proposed: a one by one increased number of filter coefficients is quantized where the augmented frequency domain error is re-minimized in each case. This new approach achieves a smaller frequency domain error than rounding of the optimal non-quantized coefficients. It is also successfully applied to the design of Lth band filters with minimum frequency domain error. FIR filters up to a filter order 1500 are designed. Keywords— Design of high order FIR filters, Constrained Least Squares (CLS) design, finite wordlength coefficients
I. Introduction Conventional design methods of FIR filters with linear phase minimize a norm of the frequency domain error E(Ω) = H(Ω) − D(Ω)
(1)
on a prescribed approximation domain B, being a subset of [0, π]. H(Ω) is the real-valued amplitude response and D(Ω) a given real-valued target function. (In this contribution no weighting function is considered.) Two error norms are commonly used: the L∞ -error norm (maximum error) ||E(Ω)||∞ = max |E(Ω)| Ω∈B
and the L2 -error norm (mean squared error) Z 1 2 ||E(Ω)||2 = |E(Ω)|2 dΩ. 2π Ω∈B
(2)
(3)
The coefficients of a filter with minimum L∞ -error norm are computed by the well known MPR-program [16]. The filter coefficients for a minimum L2 -error norm are obtained by solving a system of linear equations [3]. Both algorithms yield infinite wordlength coefficients. However the realization of a digital filter with finite expenditure allows inevitably only filter coefficients of finite precision. If a digital signal processor with floating point arithmetic is used, the effects of a finite coefficient wordlength w can (nearly always) be neglected. However for a fast and cheap implementation of a digital filter (e.g. as ASIC) a fixed point arithmetic has to be used (as assumed in the following) and the effects of finite coefficient wordlength must be considered. This is especially true for very long filters with applications to sharp cut-off or sample rate alteration [24]. The author is with SIEMENS ICM, Munich. email: [email protected]. This work was performed during his time as research assistant with Digital Signal Processing Group at the Ruhr-Universit¨ at Bochum, Germany and was supported by Deutsche Forschungsgemeinschaft under contract GO 849/1-1.
The most widespread method of finding finite wordlength coefficients is the direct quantization (DQ) method in which the infinite wordlength coefficients, obtained by one of the aforementioned filter design methods [16], [19], are rounded to yield quantized coefficients. As a result the respective error norm increases, possibly beyond the maximally allowed value. To achieve a smaller error norm the coefficient wordlength and/or the filter order must be increased [9], leading to an additional expenditure. As the finite wordlength coefficients found by the direct quantization method are not optimal in sense of the error norm (since this method does not take into account the coefficient wordlength), finite wordlength filters of original order and coefficient wordlength may exist with smaller error norm. Several contributions deal with the problem of finding the optimal finite wordlength coefficients of FIR filters with minimum L∞ -error norm by using mixed integer linear programming [10], [14], [15], simulated annealing [2] and optimal or local search methods [4], [11], [25]. All these methods have nothing in common with the conventional infinite wordlength FIR filter design method [16] and are very time-consuming or, even worse, do not converge for high filter orders. A more promising approach, where the coefficients are quantized one by one and the L2 -error norm is reminimized each time, is presented in [13]. An application of this approach to the design of FIR filters with minimum L∞ -error norm requires an algorithm for re-minimization of the L∞ -error norm. A modification of the MPR algorithm [16] can not be recommended because the MPR algorithm can only influence as many extrema as free variables (number of non-quantized filter coefficients) exist, and the L∞ error norm increases significantly with a decreasing number of extrema [6], [5]. This contribution focuses on the design of linear-phase FIR filters of high order and minimum L∞ -error norm (2). A filter with minimum L∞ -error norm is further on referred to as optimal. Instead of a modified MPR algorithm a modified Constrained Least Square (CLS)-algorithm [1] is proposed for re-minimization. The filter coefficients are assumed to be identical with the impulse response h(k) as it is the case for a direct form implementation. The coefficients are restricted to take only discrete values on the quantization grid [−1 + q, −1 + 2q, . . . , 1 − q, 1] with q = 2−w+1 . First a procedure similar to the method in [13] is outlined, where the filter coefficients are iteratively quantized. Then the modified CLS-algorithm is presented which is used for the re-minimization of the L∞ -error norm after the quantization of one somehow selected coefficient. The
2
EURASIP SIGNAL PROCESSING
design examples show the power of this new approach.
error: min |∆h(k)| = min |hQ (k) − h(k)|.
II. Iterative Quantization Procedure The proposed Iterative Quantization (IQ) consists of 3 steps. Steps 2 and 3 are repeated iteratively. In each iteration one additional coefficient is quantized. So the set of quantized coefficients increases from iteration to iteration until all coefficients are quantized. The IQ procedure is outlined subsequently (Fig. 1): 1. Initialization: The L∞ -error norm is minimized with the MPR [16] yielding the non-quantized coefficient set {h0 (k)}. The set {h0Q (kf )} of quantized coefficients and the set {kf } of indices of quantized coefficients are void. 2. Selection and Quantization of an Additional Coefficient: Out of the actual set {hi−1 (k)} of non-quantized coefficients one coefficient h(kS ) is selected, rounded to the next quantization step and assigned to the next set {hiQ (kf )} of quantized (fixed) coefficients. The index kS is shifted from the set {k} of indices of non-quantized coefficients to the set {kf } of indices of quantized coefficients. Hence, from one iteration (i−1) to the next i, the set {hiQ (kf )} increases by one coefficient, while the set {hi (k)} diminishes by one coefficient. 3. Re-Minimization of L∞ -Error Norm: The L∞ -error norm increased due to the preceding coefficient quantization is re-minimized. Note that only the set {hi (k)} can be used for re-minimization. The iterative procedure ends as soon as the set {hi (k)} is void.
∀k
III. Re-Minimization of L∞ -Error Norm For the re-minimization of the L∞ -error norm considering the set {hiQ (kf )} of fixed filter coefficients a modified CLS-algorithm is used. The CLS-algorithm [1] minimizes the L2 -error norm (3) while the L∞ -error norm (2) is constrained not to exceed prescribed bounds. To use the CLSalgorithm for the minimization of the L∞ -error norm these prescribed bounds must be chosen as tighten as possible. The combination of the Iterative Quantization with the modified CLS-algorithm for the re-minimization of the L∞ error norm is subsequently called IQCLS. First the original CLS-algorithm [1] is described and then its modification. A. CLS-Algorithm for the Design of FIR Filters The amplitude response H(Ω) of a zero-phase FIR filter of even order n and even symmetry1 is H(Ω) = h(0) +
h
= [1, 2 cos(Ω), . . . , cos( n2 Ω)], T
= [h(0), . . . , h( n2 )] .
(7)
(8)
Since E(Ω) is real and a scalar:
|E(Ω)|2 = E 2 (Ω) = E(Ω)T E(Ω).
Quantization of h(kS) hQ(kS) {hQi(kf)}
i=i+1
Re-Minimization of L -error {hi(k)}
where
Fig. 1. Flow chart of the Iterative Quantization (IQ) Procedure
There exists a large number of possible selection criteria for step 2. The most successful selection strategy [6], [5] chooses the coefficient with the smallest absolute rounding
(9)
Introducing (9) with (8) in (3) leads to an alternative expression of the L2 -error norm [3] Z 1 ||E(Ω)||22 = hT Ah − 2hT b + D2 (Ω)dΩ, (10) 2π Ω∈B
(6)
(5)
With (5) the error E(Ω) of (1) is
E(Ω) = aT h − D(Ω).
2 cos(kΩ)h(k) = aT h
with2
i=0
Selection of one h(kS) out of {hi-1(k)}: kS {kf}, {k}={k} kS
n/2 X k=1
i=1
(4)
The IQ method is suboptimal but leads with a great probability to a smaller L∞ -error norm than the direct quantization of the coefficients from the MPR-program (DQMPR), because the L∞ -error norm is constantly reminimized considering the already quantized coefficients.
aT Initialization: Minimization of L - error {h0(k)} with { kf},{hQ0(kf)} void
∀k
A
=
b
=
Z 1 aaT dΩ, 2π Ω∈B Z 1 aD(Ω)dΩ. 2π Ω∈B
(11) (12)
1 Equivalent equations can easily be found for all other types of linear phase FIR-filters [19]. 2 Boldface capital letters are used for matrices, boldface lower case letters for vectors.
EVANGELISTA: DESIGN OF DIGITAL FIR FILTERS
3
The minimization of (10) with respect to h leads to the starting solution of the CLS-algorithm: hstart = A
−1
b.
(13)
If all local extrema of the amplitude response of hstart do not exceed the prescribed upper and lower bounds Rlo (Ω) and Rup (Ω) on B, no iteration is necessary and hstart is the final solution of the CLS-algorithm. Otherwise the amplitude response is set to the upper or lower bound, respectively, at all extremal frequencies Ων , where H(Ω) violates these bounds: ½ Rup (Ων ), if H(Ων ) > Rup (Ων ), ! H(Ων ) = (14) Rlo (Ων ), if H(Ων ) < Rlo (Ων ), or in matrix notation Bh = r
(15)
with
.. . T B= a (Ων ) , .. .
.. .
r= Rup,lo (Ων ) . .. .
(16)
Hence, in this case the extended L2 -error norm, the Lagrange-function [8] ELag = ||E(Ω)||22 − mT (Bh − r)
(17)
with the Lagrange-multipliers m, has to be minimized. The minimization of the Lagrange-function (17) in conjunction with (10) with respect to h leads to Ah − 0.5B T m = b.
(18)
Eqs. (15) and (18) commonly represent a set of linear equations with the solution [1]: m = 2(BA−1 B T )−1 (r − BA−1 b) h = A−1 (b + 0.5B T m).
(19) (20)
The amplitude response of h is calculated and if its extrema are beyond the bounds at the (new) extremal frequencies Ων , the procedure starting at eq. (14) is repeated. The CLS-algorithm stops if all extrema of the amplitude response are within the bounds (convergence) or the set {Ων } of extremal frequencies does not change anymore (no convergence, because the bounds are to tight). B. Modification of the CLS-Algorithm To consider fixed coefficients hQ (kf ) in the CLSalgorithm the amplitude response (5) is expressed depending of a vector hf containing the fixed coefficients hQ (kf ) and a vector hu containing all non-quantized coefficients to be used for re-minimization: T H(Ω) = aT u hu + af hf .
(21)
Replacing H(Ω) with (21) in (1) yields E(Ω) = aT u hu − DMod (Ω),
(22)
DMod (Ω) = D(Ω) − aT f hf .
(23)
where As the expressions of the modified error (22) and the original error (1) and (5) are the same replacing a with au , h with hu and D(Ω) with DMod (Ω), the modified and original CLS-algorithm are also the same by using the same replacements. So A (11), b (12) and B (16) must be replaced with Z 1 (24) au aT AMod = u dΩ, 2π Ω∈B Z 1 au DMod (Ω)dΩ, (25) bMod = 2π Ω∈B .. T. BMod = (26) au (Ων ) . .. . The matrix AMod to be inverted in the modified CLSalgorithm is of lower order than the matrix A of the original CLS-algorithm. So the consideration of the fixed coefficients hf leads to a saving in computational expenditure. The upper and lower bounds are introduced as follows Rup (Ω) = Rlo (Ω) =
D(Ω) + δ, D(Ω) − δ,
(27) (28)
with δ as input parameter of the modified CLS-algorithm. δ must be iteratively decreased for the re-minimization of the L∞ -error norm (step 3 of the IQ) until the modified CLS-algorithm does not converge anymore. The smallest value of δ with convergence of the modified CLS-algorithm is the minimal L∞ -error norm. A major advantage of the IQCLS algorithm is that is based on a conventional design method (the CLS-algorithm [1]): All experience, especially with the design of filters of high order, can be applied. IV. Examples This section focuses on the design of lowpass filters, because lowpass filters are widely used and similar results are also obtained for all other types of frequency selective filters. A first target function D(Ω) is ½ L, for 0 ≤ Ω ≤ ΩP , D(Ω) = (29) 0, for ΩS ≤ Ω ≤ π, where the factor L in D(Ω) is chosen in such a way that the quantization grid [−1 + q, −1 + 2q, . . . , 1 − q, 1] for the coefficients hQ (kf ) is optimally exploited. In Fig. 2 the magnitude resonse |H(Ω)/L| of a lowpass filter designed with IQCLS is depicted. Although already 7 of 11 coefficients are quantized with 8 bit, |H(Ω)/L| shows nearly equi-ripple behaviour.
4
EURASIP SIGNAL PROCESSING
3) outperforms the DQMPR approach, if the difference between the minimax errors of the DQMPR and MPR is large (see figures 4, 5, 6).
Magnitude / dB
0.2
0.1
0
−0.1
−10 −0.2
0
0.05
0.1
0.15
0.2
Ω/π
0.25
0.3
0.35
0.4
−15
L −error norm / dB
−35
−40
−45
−50 0.5
0.6
0.7
0.8
0.9
1
Ω/π
1.1
1.2
1.3
1.4
1.5
−20
−25
∞
Magnitude / dB
−30
−30
Fig. 2. 20 log10 |H(Ω)/L| with IQCLS for n = 20, ΩP = 0.4π, ΩS = 0.6π, 7 of 11 coefficients quantized with w = 8: passband (above), stopband (below)
−35 50
100
150
200
Filter order n
A. Design of FIR Lowpass Filters FIR lowpass filters were designed for different specifications (29) and wordlengths w using the IQCLS algorithm. The following observations can be reported: −13
n=50
• For a given wordlength w and specification D(Ω) a filter order nmax exists beyond which the L∞ -error norm does not improve anymore for both methods, IQCLS and DQPMR (see figure 5: nmax ≈ 250 for w = 8, figure 6: nmax ≈ 180 for w = 6). This result is consisting with [9].
n=60
−15
−10
−16
n=70 −20
−17
w=8 −30
n=80
−18
−40 −19
L∞−error norm / dB
∞
L −error norm / dB
−14
Fig. 4. Log. L∞ -error norm with IQCLS (◦), DQMPR (∗) and MPR (+) for ΩP = 0.192π, ΩS = 0.208π, L = 5, w = 8
n=90
−20
n=100 −21 0
5
10
15
20
25
30
35
40
45
50
Number of quantized coefficients
−60
w=16
−70 −80 −90
−100
Fig. 3. Remaining log. L∞ -error norm after IQCLS re-minimization for ΩP = 0.192π, ΩS = 0.208π, L = 5, w = 8
−110 −120
For a small number of quantized coefficients the remaining L∞ -error norm increases scarcely after re-minimization in step 3 of the IQCLS algorithm (s. figure 3). This is due to the fact that for IQCLS re-minimization enough nonquantized coefficients are left (e.g. for n = 90, 100 in figure 3). In the opposite case a smaller final L∞ -error norm can sometimes be achieved by quantizing all remaining nonquantized coefficients directly. The number of coefficients to be quantized in one step can be deduced of the L∞ -error norm of figure 3 (e.g. for n = 50, 70, 80 as indicated by the dashed lines). • The overall increase of the L∞ -error norm during the IQCLS is smaller for larger coefficient wordlength w. • As expected the IQCLS algorithm (with adequate direct quantization of the remaining coefficients as in figure
w=12
−50
•
100
200
300
400
500
600
700
800
900
Filter order n
Fig. 5. Log. L∞ -error norm with IQCLS (◦), DQMPR (∗) and MPR (+) for ΩP = 0.192π, ΩS = 0.208π, L = 5
FIR filters were successfully designed using the IQCLS algorithm (implemented in Matlab) up to a filter order of n = 1500. Beyond this upper bound for n the Choleskydecomposition [7], [8] used to invert AMod does not work properly. •
The IQCLS method needs considerably more computation time than the DQMPR (e.g. 750 re-minimizations of L∞ error are required for n = 1500 in IQCLS while in DQMPR there is only one minimization of L∞ -error).
EVANGELISTA: DESIGN OF DIGITAL FIR FILTERS
5
−10
−29.5
−30
L −error norm / dB
−30.5
−20
−25
−31
−31.5
−32
−32.5
∞
∞
L −error norm / dB
−15
−30
−33
−33.5
−34
−35 50
100
150
−34.5
200
Filter order n
The IQCLS algorithm can also be applied to the design of Lth-band (Nyquist-)filters [17] with minimum L∞ -error norm. The zero phase impulse response h(k) of Lth-band filters must hold ½ 1, for ν = 0, h(νL) = (30) 0, for ν ∈ Z{0}. Lth-band filters are advantageously used as anti-imaging filters in interpolators with upsampling factors L [24]. An appropriate target function D(Ω) for Lth-band filters is given by [22]: L, 0 ≤ Ω ≤ ΩP , 0, 2π L ν − Ωp ≤ Ω ≤
4
5
6
7
8
9
10
Fig. 7. Log. L∞ -error norm with IQCLS, (◦), CMPR (∗) and MPR π (+) for ΩP = 0.96 L , n = 40L, non-quantized coefficients
V. Conclusion and Outlook
B. Design of Lth-Band Filters
½
3
Upsampling factor L
Fig. 6. Log. L∞ -error norm with IQCLS (◦), DQMPR (∗) and MPR (+) for ΩP = 0.192π, ΩS = 0.208π, L = 5, w = 6
D(Ω) =
2
2π L ν
+ Ωp ,
(31)
with ν = 1, . . . , L−1. Obviously D(Ω) has L−1 stopbands. The task of designing Lth-band filters with a minimum L∞ -error norm is addressed by [2], [12], [15], [20], [21] using methods with high computational expenditure. The use of the IQCLS algorithm makes the design very fast: The coefficients given by (30) are all treated as fixed coefficients in step 2 and the re-minimization of the L∞ -error norm by the modified CLS algorithm is done once in step 3 yielding the optimal, non-quantized coefficients. Of course, the design procedure can be continued in order to have all coefficients quantized, as treated above. The L∞ -error norm of Lth-band filters for different upsampling factors L and design methods is depicted in figure 7. The coefficients h(k) resulting from the MPR-program with D(Ω) according to (31) violate (30). Hence the result is not an Lth band filter. The corresponding L∞ -error norm is only depicted as lower bound of the achievable L∞ error norm. The forced correction of the coefficients h(k) from the MPR to hold (30) is done by the CMPR. As can clearly be seen the IQCLS algorithm leads to a considerably smaller L∞ -error norm than the CMPR [6].
The design of linear phase FIR filters with finite wordlength coefficients and minimum L∞ -error norm using an Iterative Quantization approach combined with a modified CLS algorithm (IQCLS) has been presented. The IQCLS algorithm leads to a smaller L∞ -error norm than the DQMPR (direct quantization of the coefficients h(k) resulting from the MPR program) for the case, where a great difference between the L∞ -error norms of MPR and DQMPR exists. The IQCLS has successfully been applied to filters up to a filter order of n = 1500. It is also very suitable for the design of Lth-band filters with minimum L∞ -error norm. For a further reduction of the L∞ -error norm the direct quantization of one or more remaining coefficients at the end of the IQCLS (see sect. II) might be replaced by applying local search algorithms [11]. The IQCLS can also be easily applied to the design of filters with different coefficient wordlengths. This leads to simple and fast multiplications for those coefficients with very small wordlengths. Especially the coefficients quantized in the first iteration steps of the IQCLS require only a small wordlength [6] (see also figure 3). The design of filter with solely power-of-two coefficients or with coefficients in CSD-code representation [18] is possible, too. With a small modification, the IQCLS can also be used for the design of finite wordlength linear phase FIR filters and Lth-band filters with minimum L2 -error norm. Hence, with the IQCLS a very flexible and powerful filter design tool is given. Because the selection criterion (4) used in the second step of the IQ is not proven to be optimal, further research could deal with finding other more appropriate selection criteria. Acknowledgement The author wishes to express his appreciation to Prof. Dr.-Ing. H. G. G¨ockler and the unknown reviewers for their careful review of the manuscript and valuable suggestions.
6
EURASIP SIGNAL PROCESSING
References [1] [2]
[3] [4]
[5] [6] [7] [8] [9] [10]
[11]
[12]
[13] [14] [15] [16]
[17] [18] [19] [20] [21] [22] [23]
[24] [25]
J. W. Adams, FIR Digital Filters with Least Squares Stopbands Subject to Peak-Gain Constraints, IEEE Transactions on Circuits and Systems, 39:376-388, April 1991. N. Benvenuto, M. Marchesi, and Aurelio Uncini, Applications of Simulated Annealing for the Design of Special Digitals Filters, IEEE Transactions on Signal Processing, 40:323-332, February 1992. C. S. Burrus, A. W. Soewito, R. A. Gopinath, Least Squared Error FIR Filter Design with Transition Bands, IEEE Transactions on Signal Processing, 40:1327-1340, June 1992. C.-L. Chen and A. N. Willson Jr., A Trellis Search Algorithm for the Design of FIR Filters with Signed-Powers-of-Two Coefficients, IEEE Transactions on Circuits and Systems II, 46:2939, January 1999. M. Erdogan, Vergleichende Untersuchung zum Entwurf linearphasiger FIR-Filter mit quantisierten Koeffizienten, Diplomarbeit, Ruhr-Universit¨ at Bochum, 1999. G. Evangelista, Zum Entwurf digitaler Systeme zur asynchronen Abtastratenumsetzung, Ph.D. Thesis, Bochum, 2000. F. R. Gantmacher, Matrizentheorie, Berlin: Springer, 1986. G. H¨ ammmerlin and K.-H. Hoffman, Numerische Mathematik, Berlin: Springer, 1991. ¨ U. Heute, Uber Realisierungsprobleme bei nicht-rekursiven Digitalfiltern, Ph.D. Thesis, Erlangen, 1975. D. M. Kodek, Design of Optimal Finite Wordlength FIR Digital Filters using Integer Programming Techniques, IEEE Transactions on Acoustics, Speech, Signal Processing, 28:304-308, June 1980. D. M. Kodek and K. Steiglitz, Comparison of Optimal and Local Search Methods for Designing Finite Wordlength FIR Digital Filters, IEEE Transactions on Circuits and Systems, 28:28-32, January 1981. J. K. Liang, R. J. P. DeFigueiredo, and F. C. Lu, Design of Optimal Nyquist, Partial Response, N th Band, and Nonuniform Tap Spacing FIR Filters using Linear Programming Techniques, IEEE Transactions on Circuits and Systems, 32:386392, April 1985. Y. C. Lim and S. R. Parker, A Discrete Coefficient FIR Digital Filter Design Based Upon an LMS Criteria, Int. Symposium Circuits and Systems, Rome, Italy, 3:796-799, May 1982. Y. C. Lim and S. R. Parker, FIR Filter Design Over a Discrete Powers-of-Two Coefficient Space, IEEE Transactions on Acoustics, Speech, Signal Processing, 31:583-591, June 1983. Y. C. Lim and B. Liu, Design of Cascade FIR Filters with Discrete Valued Coefficients, IEEE Transactions on Acoustics, Speech, Signal Processing, 36:1735-1739, November 1988. J. H. McClellan, T. W. Parks and L. R. Rabiner, A Computer Program for Designing Optimum FIR Linear Phase Digital Filters, IEEE Transactions on Audio Electroacoustics, 21:506-526, 1973. F. Mintzer, On Half-band, Third-Band, and N th Band FIRFilters and Their Design, IEEE Trans. on Acoustics, Speech, Signal Processing, ASSP-30:734-738, October 1982. S. K. Mitra and J. F. Kaiser, Handbook for Digital Signal Processing, p. 206, J. Wiley & Sons: New York, 1993. T. W. Parks and C. S. Burrus, Digital Filter Design, New York: John Wiley & Sons, 1987. H. Samueli, On the Design of Optimal Equiripple FIR Digital Filters for Data Transmission Applications, IEEE Transactions on Circuits and Systems, 35:1542-1546, December 1988. T. Saram¨ aki and Y. Neuvo, A Class of FIR Nyquist (N thBand) Filters with Zero Intersymbol Interference, IEEE Transactions on Circuits and Systems, 34:1182-1190, October 1987. R. W. Schafer and L. R. Rabiner, A Digital Signal Processing Approach to Interpolation, Proceedings of the IEEE, 61:692702, June 1973. I. W. Selesnick, M. Lang and C. S. Burrus, Constrained Least Square Design of FIR Filters without Specified Transition Bands, IEEE Transactions on Signal Processing, 44:1879-1892, August 1996. P. P. Vaidyanathan, Multirate Systems and Filter Banks, Englewood Cliffs: Prentice Hall, 1993. Q. Zhao and Y. Tadakoro, A Simple Design of FIR Filters with Powers-of-Two Coefficients, IEEE Transactions on Circuits and Systems, 35:566-570, May 1988.
