* Corresponding author.
E-mail address: sgchen@cc.nctu.edu.tw (Sau-Gee Chen).
Design of "nite-word-length FIR "lters with least-squares error
Yung-An Kao, Sau-Gee Chen*
Department of Electronics Engineering and Institute of Electronics, National Chiao Tung University, 1001 Ta Hsueh Road, Hsinchu, Taiwan
Received 25 May 1998; received in revised form 5 January 2000
Abstract
This paper proposes a new algorithm for designing "nite word length linear-phase FIR "lters. The new algorithm produces "nite-precision least-squares error (LSE) solutions with much reduced search time than the brute-force full search algorithm. It is di!erent from the full search algorithm that tries all possible combinations directly. The new algorithm utilizes geometric properties of a hyper-space to pinpoint potential solutions in a much more restricted way. Accordingly, a much smaller search space is generated. 2000 Elsevier Science B.V. All rights reserved.
Zusammenfassung
In dieser Arbeit wird ein neuer Entwurfsalgorithmus fuKr linearphasige FIR-Filter bei endlicher WortlaKnge vorges-chlagen. Der neue Algorithmus liefert LSE (kleinstes Fehlerquadrat)-LoKsungen mit endlicher Genauigkeit bei sehr verkleinerter Suchdauer gegenuKber der vollstaKndigen Suche. Er ist verschieden von einer vollstaKndigen Suche, die alle moKglichen Kombinationen direkt ausprobiert. Der neue Algorithmus nutzt geometrische Eigenschaften eines Hyper-raumes aus, um potentielle LoKsungen in einer eingeschraKnkten Weise festzulegen. Dadurch wird ein viel kleinerer Suchraum erzeugt. 2000 Elsevier Science B.V. All rights reserved.
Re2sume2
Nous proposons dans cet article un algorithme nouveau pour la conception de "ltres FIR a` phase lineHaire en preHcision "nie. Cet algorithme produit des solutions aux monidres carreHs (LES) avec un temps de recherche bien plus reHduit que l'approche de recherche exhaustive. II est di!eHrent de l'alogorithme de recherche exhaustive qui essaye directement toutes les combinaisons possibles. Cet algorithme utilise les proprieHteHs geHomeHtriques d'un hyperspace pour mettre en eHvidence less solutions potentielles d'une manie`re beaucoup plus restrictive. De ce fait un espace de recherche beaucoup plus petit est geHneHreH. 2000 Elsevier Science B.V. All rights reserved.
1. Introduction
In practice, "lters are realized by "xed-point arithmetic. In designing "nite word length or
powers-of-two linear-phase FIR "lters there are many algorithms based on integer programming [5] and the modi"ed integer programming rithms [4,6}8,10,11,13]. Solutions of these algo-rithms are found by searching the regions con"ned by some linear constraints subject to minimizing objective functions. The computation load of the linear/integer programming approach [5] is very
0165-1684/00/$ - see front matter 2000 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 0 0 ) 0 0 0 5 8 - X
Nomenclature
u frequency in radian
D(u) desired frequency response
uN passband cuto! frequency uQ stopband cuto! frequency
h(n) "lter coe$cient
h optimal coe$cient vector
D system matrix for solving h
C vector for solving h, i.e., h"D\C h Rounded vector of h
E(h) square error function due to coe$cient vector h
h optimal "nite-precision coe$cient vector producing the least-squares error
h "nite-precision version of the coe$cient vector h
h(n) "nite-precision version of the coe$cient h(n) h (n) lower bound of h(n)
h(n) upper bound of h(n)
S(n) number of all the candidate "nite-precision coe$cients of h(n) for the full search algo-rithm
¸(n) number of all the candidate "nite-precision coe$cients of h(n) for the new optimization algorithm
heavy, and it is intended for minimizing the min}max error norm. In implementation, all these algorithms need to sample the target "lter spectrum for testing constraints, instead of ideally testing the whole continuous frequency band. This results in computation penalty, as well as error. The algo-rithms in [4,10] provided fast search algoalgo-rithms to reduce computation time. Some of the local search algorithms [6,13] reduce search time, at the expense of performance. There are the e!ective but computa-tion-intensive simulated annealing technique [1,2]. Simulated annealing methods require very intensive computation. The near least-squares error ap-proaches [9,12] reduce computation time consider-ably, but only get the suboptimal solutions.
In summary, the existing algorithms either pro-duce optimal results at the cost of intensive compu-tation, or produce suboptimal results at a much reduced computation load. In this paper, we will propose a new LSE optimization algorithm for "nite word length "lters. For each coe$cient, the new algorithm utilizes the geometric projection property of a hyper-space to locate potentially dis-crete solutions, subject to the LSE constraint. From these possible solutions, an e$cient tree path search method is introduced to pinpoint the "nal optimal LSE solution. Doing this way, a much smaller search space than that of brute-force search algorithm is generated, and accordingly a much reduced search time.
2. The new algorithm
The new algorithm starts with the optimal in"-nite-precision LSE solution [3] to the given ideal response D(u), where D(u)"1, for 0)u)uN,
D(u)"0, for u)u)n. Without loss of generality,
we consider an N-tap, symmetric, zero-phase, odd-length "lter, with the frequency response, H(eS)"
h(0)#2 ,\
L h(n) cos (nu), h(n)"h(!n). The
op-timal LSE solution h"[h(0)h(1)2h((N!1)/2)]2 can be solved as h"D\C [3], by setting the gradient of the square-error cost function
E(h)"p uN
SN [D(u)!H(eS)] du #s uQ L SQ [H(eS)] duto zero, where C and D are vector and matrix, respectively, depending uponuN, uQ, p and s.
De"ne h"round(h2@)2\@, with b#1 the num-ber of bits, we can get a square error k"E(h). Note that k is very close to the least-squares error
E(h) and the error surface E(h)"k encloses h. As
will be shown later in the simulations h is a good initial value for locating the optimal discrete
h and occasionally h"h. Therefore, we can
"nd some discrete coe$cient vectors whose square errors are smaller than k if they are inside the error
surface E(h)"k. On the other hand, h is the opti-mal discrete solution when there is no discrete coe$cient vector inside the error surface E(h)"k. The design problem then reduces to: how do we locate these discrete points which are inside the error surface E(h)"k in an e$cient way? To solve the design problem, we will iteratively use a projec-tion algorithm in "nding potential discrete coe$-cients, in combination with an e$cient tree-path search algorithm. The projection algorithm utilizes the geometric properties of an LSE surface. Before introducing the new algorithm, we "rst introduce the projection algorithm.
2.1. The projection algorithm
Given a hyper-ellipse described by E(h)"k, if there exits a discrete coe$cient vector enclosed by the hyper-ellipse, then the coe$cient vector will produce a square error smaller than k. From geo-metric point of view, to "nd all the potential "nite-precision solutions of a particular coe$cient h(m), one can project the hyper-ellipse onto the h(m) axis. This results in a line segment enclosed by
h(m)"h (m) and h(m)"h(m), h (m)(h(m). All
the discrete h(m) points within the line segment potentially lead to a smaller square error than k. Obviously, the projection process is done by locat-ing two surfaces tangent to the hyper-ellipse. Geo-metrically, the projection is required to be tangent to the hyper-ellipse, and parallel to all h(n) axes,
n"0,2, (N!1)/2 and nOm, but perpendicular
to the h(m) axis. Hence, the two tangent points must satisfy the condition that *E(h)/*h(n)"0,
n"0,2, (N!1)/2, nOm. The condition results
in a set of (N!1)/2 linear equations. From these equations, the coe$cients h(0),2, h((N!1)/2) excluding h(m) can be solved in terms of h(m). That is, they can be solved as hK"(DK)\CK, where
hK is the coe$cient vector excluding h(m), DK is
the system matrix of the set of linear equations, and
CK is a vector whose elements are composed of
linear combinations of h(m) and constants. Since the tangent points are on the hyper-ellipse, we can substitute all h(n)'s, which are all linear functions of
h(m), n"0,2, (N!1)/2, nOm, into the quadratic
hyper-ellipse function E(h)"k. As a result, we have a quadratic equation of h(m) whose roots are h (m)
and h(m), which are the end points of the projected line segment of the hyper-ellipse. In between these two points there are S(m) discrete values of h(m).
2.2. The new xnite-precision LSE algorithm based on the projection algorithm
Assume that the "nite-precision solutions
h(0),2, h(m!1) for coe$cients h(0),2, h(m!1) have been temporarily found and "xed in
a manner as described in the following treatment for h(m) of h(m) similar to the projection method introduced before, then all the potential "nite-pre-cision LSE solutions h(m)'s for h(m) can be found by setting the gradient of E(h) to zero as
*E(h)
*h(n)"0, n"m#1,2, (N!1)/2,
which results in a set of (N!1)/2!m linear equa-tions. From these equations, the coe$cients
h(m#1),2, h((N!1)/2) can be solved in terms of h(m). By plugging these solutions into equation E(h)"k, one can solve two real roots h (m) and h(m) of h(m), h (m)(h(m). In between these two
points there are ¸(m) discrete values of h(m). By combining the projection algorithm iteratively with an e$cient tree search scheme, one has the follow-ing optimization algorithm.
2.2.1. The optimization process of the new algorithm Step 1. Solve the in"nite-precision LSE solution
h.
Step 2. Get h by directly rounding h. Let m"0, E "k"E(h) and let the optimal
discrete solution h"h.
Step 3. Find all the ¸(m) potential discrete values h(m) of h(m) using the projection
algo-rithm. Reset the index j(m) (of the candidate discrete values) of h(m) to j(m)"0. Note that all the coe$cients h(0),2, h(m!1) here have been replaced with some discrete values in the cost function E(h).
Step 4. Let j(m)"j(m)#1. Replace h(m) with its j(m)th discrete value in E(h), where h is the
coe$cient vector consisting of h(0),2,
Table 1
Speed and LSE comparisons between new algorithm and the Shyu and Lin algorithm [12], and the full search algorithm Filter
length
Square error Computation time (s) No. of all possible solutions
for the full search algorithm New algorithm Algorithm of [12] Due to h New algorithm Algorithm of [12]
*19 9.6157e!4 9.6157e!4 9.6157e!4 2.100e!1 1.800e!1 3.2000e#1
21 4.3243e!4 4.3243e!4 4.3244e!4 3.200e!1 2.000e!1 4.6080e#3
23 4.3077e!4 4.3077e!4 4.3086e!4 1.630e#0 2.300e!1 3.1104e#5
25 2.0633e!4 2.0634e!4 2.0637e!4 1.560e#0 2.700e!1 3.4992e#5
*27 1.8500e!4 1.8500e!4 1.8500e!4 5.600e!1 3.100e!1 1.8662e#5
29 1.0632e!4 1.0632e!4 1.0636e!4 1.860e#0 3.500e!1 3.7791e#7
31 7.7723e!5 7.7723e!5 7.7734e!5 8.900e!1 4.000e!1 1.0078e#7
*33 5.6581e!5 5.6581e!5 5.6581e!5 1.080e#0 4.900e!1 1.7916e#8
35 3.3187e!5 3.3187e!5 3.3221e!5 2.110e#0 5.700e!1 3.2249e#9
37 2.9730e!5 2.9730e!5 2.9806e!5 7.900e#0 6.400e!1 5.6435e#11
39 1.5046e!5 1.5103e!5 1.5127e!5 2.847e#1 7.200e!1 2.4079e#13
*41 1.4977e!5 1.4977e!5 1.4977e!5 5.170e#0 8.100e!1 9.2096e#12
43 7.3998e!6 7.4425e!6 7.4665e!6 2.600e#1 9.100e!1 8.2591e#14
45 7.1164e!6 7.1331e!6 7.2659e!6 2.274e#2 1.020e#0 2.8400e#18
47 3.9814e!6 3.9957e!6 4.0026e!6 2.529e#1 1.140e#0 2.5142e#17
49 3.2726e!6 3.3374e!6 3.3374e!6 3.657e#1 1.260e#0 5.8078e#19
51 2.2629e!6 2.2955e!6 2.3479e!6 2.301e#2 1.370e#0 1.0061e#23
discrete values of h(m#1),2, h((N!1)/2) remain to be determined.
Step 5. Cases:
(i) h(m) is not the last coe$cient and at least
one of the discrete values of h(m) has not been tested (that is, j(m))¸(m)), let
m"m#1 and go to Step 3.
(ii) h(m) is the last coe$cient and at least one of
the discrete values of h(m) has not been tested (that is, j(m))¸(m)), then the bottom level is reached and a complete discrete vector h is obtained, do the operations:
h"h and E "E(h) if E(h)(E ,
go to Step 4.
(iii) Here, all the discrete values of h(m) have been tested (that is, j(m)'¸(m)). Let
m"m!1, go to Step 4 if m*0 (regardless
of whether h(m) is the last coe$cient or not), otherwise go to step 6.
Step 6. All the h's have been searched and the LSE
solution is obtained, end the optimization process.
According to simulations, most of search paths did not go to the bottom coe$cient level, because
in most cases the projection algorithm produces null discrete solutions in higher levels. This prop-erty greatly reduces the optimization time.
3. Simulations
A low-pass "lter design problem is simulated. All the simulations were performed on UltraSPARC using MATLAB 5.1. Here, the "lter length N is varied from 19 to 51 (where N is an odd num-ber), uN"0.4p, uQ"0.5p, p"0.5, s"0.5, and wordlength"12 bits. The detail simulation data is summarized in Table 1, where the mark &*' indicates the cases when h"h. In the table, we only list the numbers of all possible solutions for the full search algorithm, because the computation times of full search algorithm greatly increase with N and far exceed those required by non-full search algo-rithms. The full search algorithm is also based on the projection algorithm de"ned in subsection 2.1 of Section 2. Speci"cally, there are S(0)S(1)2
S((N!1)/2!1)S((N!1)/2) coe$cient vectors to
be simulated. The number of combinations in-creases exponentially with "lter length.
Fig. 1. Frequency response comparison, "lter length"39, word length"12 bits.
To compare the new optimal algorithm with the existing e$cient (however, non-optimal) rithms, we simulated the fast but suboptimal algo-rithm by Shyu and Lin [12] (which we consider the most e$cient algorithm in the literature), using the same design example. Parameter ¸ in [12] is set to 3. The frequency responses of new algorithm and the algorithm of [12] are shown in Fig. 1 for
N"39. As shown, for smaller N, Shyu and Lin's
algorithm can obtain the same optimal results as those of the new algorithm in most cases, within shorter time duration than those of the new rithm. However, for larger N, Shyu and Lin's algo-rithm fail to locate the optimal solutions. Also, notice that, for o!-line and "xed-coe$cient ap-plications, "lter design time is generally not an issue as long as one can "nd the optimal solution within an allowable amount of time. This argument puts the new algorithm in a more appealing position than the highly cost-e!ective (but suboptimal) algo-rithm of [12].
Table 1 also shows the square errors due to h. As shown, h's are good initial values for locating the optimal h's, that give square errors close to the LSE's produced by h's. In some cases h is equal to h. In this situation, the new algorithm can solve the optimal solution very quickly. As can be seen, the square errors generally reduce and computation times increase with the increasing
"l-ter length. To roughly compare the min-max ap-proach [5], we also simulated the example of [5] with the speci"cations: uN"0.4p, uQ"0.5p,
N"21, and wordlength"6 bits. In this case,
h"h and E(h)"E(h)"7.7119e!4 which
is predictably smaller than the square error
E(h } )"12.8863e!4 due to the min-max
solution h } from [5]. On the other hand, the
max error due to h is 0.1094, which is also predict-ably larger than the min-max error 0.0711 due to
h } . For other design examples, similar
com-parison results can be concluded as this one.
4. Conclusion
An e$cient "nite-precision "lter optimization al-gorithm generating LSE results is proposed. It is di!erent from the brute-force search algorithm that tries all possible combinations directly. The new algorithm utilizes geometric properties of a hyper space to pinpoint potential solutions in a much more restricted way, and accordingly a much small-er search space is gensmall-erated. The future work is to extend the algorithm to weighted LSE "lters and 2-D "lters.
References
[1] N. Benvenuto, M. Marchesi, Digital "lters design by simulated annealing, IEEE Trans. Circuits Systems 36 (March 1989) 459}460.
[2] N. Benvenuto, M. Marchesi, A. Uncini, Applications of simulated annealing for the design of special digital "lters, IEEE Trans. Signal Process. 40 (February 1992) 323}332. [3] S. Burrus, A.W. Soewito, R.A. Gopinath, Least squared error FIR "lter design with transition bands, IEEE Trans. Signal Process. 40 (6) (June 1992) 1327}1340.
[4] B. Jaumard, M. Minoux, P. Siohan, Finite precision design of "lters using a convexity property, IEEE Trans. Acoust. Speech, Signal Process. 36 (3) (March 1988) 407}411. [5] D.M. Kodek, Design of optimal "nite word length FIR
digital "lters using integer programming technique, IEEE Trans. Acoust. Speech, Signal Process. ASSP-28 (3) (June 1980) 304}308.
[6] D.M. Kodek, K. Steiglitz, Comparison of optimal and local search methods for designing "nite word length FIR digital "lters, IEEE Trans. Circuits Systems CAS-28 (1) (January 1981) 28}32.
[7] Y.C. Lim, S.R. Parker, Finite word length FIR "lter design using integer programming over a discrete coe$cient space, IEEE Trans. Acoust. Speech, Signal Process. ASSP-30 (4) (August 1982) 661}664.
[8] Y.C. Lim, S.R. Parker, FIR "lter design over a discrete powers-of-two coe$cient space, IEEE Trans. Acoust. Speech, Signal Process. ASSP-31 (3) (June 1983) 583}591. [9] Y.C. Lim, S.R. Parker, Discrete coe$cient FIR digital "lter design based upon an LMS criteria, IEEE Trans. Circuits Systems CAS-30 (October 1983) 723}739.
[10] J.P. Marques, A new design method of optimal "nite word length linear phase FIR "lters, IEEE Trans. Acoust.
Speech, Signal Process. ASSP-31 (4) (August 1983) 1032}1034.
[11] H. Samueli, An improved search algorithm for the design of multiplierless FIR "lters with powers-of-two coe$-cients, IEEE Trans. Circuits Systems 36 (7) (July 1989) 1044}1047.
[12] J.J. Shyu, Y.C. Lin, A new approach to the design of discrete coe$cient FIR digital "lters, IEEE Trans. Signal Process. 43 (January 1995) 310}314.
[13] Q. Zho, Y. Tadokoro, A simple design of FIR "lters with powers-of-two coe$cients, IEEE Trans. Circuits Systems 35 (May 1988) 566}570.
