- Email: [email protected]
Design of linear or minimum-phase fir filters by constrained chebyshev approximation
Signal Processing 5 (1983) 325-332 North-Holland
325
D E S I G N OF L I N E A R OR M I N I M U M - P H A S E FIR FILTERS BY C O N S T R A I N E D C H E B Y S H E V A P P R O X I M A T I O N F. G R E N E Z , m e m b e r E U R A S I P
Service d'Electricit~ Ggndrale, C.P. 165, UniversitdLibre de Bruxelles, av. F. Roosevelt 50, 1050 Bruxelles, Belgium Received 15 December 1982 Revised 17 February 1983
Abstract. Constrained Chebyshev approximation is a known mathematical problem, and Remes' type algorithms for finding the solution have been described. It is shown in this paper that the classical Parks-McClellan program can be converted into a general program for constrained approximation. The program can serve for the design of both linear and minimum phase FIR digital filters.
Zusammeniassung.Die Tschebyshceff-Approximation mit Nebenbedigungen ist ein bekanntes mathematicsches Problem, welches mit Hilfe von Remes-artigen Algorithmen gelSst werden kann. Es wird in diesem Beitrag gezeigt, wie das klassische Program yon Parks und McClellan in einem allgemeinen Programm fiir die Approximation mit Nebenbedingungen modifiziert werden kann. Ein solches Programm kann zur Berechnung nichtrekursiver Filter mit linearer oder mit minimaler Phase gebraucht werden. R6sum6. L'approximation de Tchebycheff avec contraintes est un probl~me math6matique connu dont on peut obtenir la solution par des algorithmes de type Remez. On montre dans cet article comment le programme classique de Parks et McClellan peut ~tre converti en un programme g6n6ral d'approximation avec contraintes. Un tel programme peut servir pour le calcul des filtres non r6cursifs ~ phase lin6aire ou h phase minimale. Keywords.Constrained approximation, FIR filter.
1. Introduction
constraints of the form:
T h e C h e b y s h e v , or m i n i m a x , a p p r o x i m a t i o n of a c o n t i n u o u s f u n c t i o n f ( x ) by p o l y n o m i a l s is a classical p r o b l e m , a n d consists of finding the polyn o m i a l p ( x ) which m i n i m i z e s the u n i f o r m n o r m of the error f u n c t i o n :
[If-pll = max x
w(x)lf(x)-p(x)l.
T h e design of a linear phase F I R digital filter can be f o r m u l a t e d as a C h e b y s h e v a p p r o x i m a t i o n p r o b l e m a n d an efficient c o m p u t a t i o n p r o g r a m is available [3]. In some applications, it can be useful to restrict the class of a p p r o x i m a t i n g p o l y n o m i a l s by a d d i n g
l(x)<~p(x)<~u(x). F o r example, the allowable d e v i a t i o n in s o m e interval can be fixed to a given value, or the a p p r o x i m a t i n g p o l y n o m i a l can be r e q u i r e d to be nonnegative. This p r o b l e m is k n o w n in the literature as the c o n s t r a i n e d or restricted a p p r o x i m a t i o n . It is recalled here how the best restricted a p p r o x i m a t i o n can b e calculated a n d it is s h o w n that, for the design of F I R digital or C C D filters, the available p r o g r a m for u n c o n s t r a i n e d approxim a t i o n can b e modified to solve the m o r e g e n e r a l p r o b l e m of c o n s t r a i n e d a p p r o x i m a t i o n .
0165-1684/83/$3.00 © 1983, Elsevier Science Publishers B.V. (North-Holland)
F. Grenez / Constrained approximation for FIR filters
326
2. C o n s t r a i n e d C h e b y s h e v a p p r o x i m a t i o n
in Xp satisfying rr(xi) = (-1)'+Lo'(xl) where ~r(x) = - 1 if x ~ X 1 u X 2 and ~ r ( x ) = + l if x
Definitions
X + a k-) X + 2 .
Let X be a compact subset of [a, b ], qh • • • ~0, a Chebyshev set on X, and M the set of generalized approximating polynomials:
p (x) = ~ ~j~0;(x). i=l
Let f(x) be the function to be approximated and u (x), l(x) two auxiliary functions (called the upper and lower bound), with:
In other words, the best constrained approximation polynomial p (x) must exhibit n + 1 extremal points, where p(x) is alternatively above and below fix), and where either the weighted error function is m a x i m u m or where the polynomial touches one bound. This theorem generalizes the classical alternation theorem for unconstrained approximation.
Algorithm [2] l(x)
x6S. l(x)<~f(x)<~u(x), The subset K of M is defined as:
K = {p ~ M: l(x) <-p(x)<~ u (x), x ~ X}. K is the set of approximating polynomials whose ranges are restricted by I and u. Then the polynomial p*(x)~ K is the best constrained Chebyshev approximation to f if
II/-p*ll
= min p~K
A R e m e s ' type algorithm for the construction of the best constrained approximation can be described as follows. Step 1. Choose n + 1 extremai points x~ < x ~ < 1 • • • < x , + l and solve the n + 1 equations: 1
1
l
~ lai o~n+l - f ( x ] ) ,
) - , - , w(-U, )
/'=1
i=1 ..... n+l for an+l. 1 1 I. The polynomial Define e 1 = la.+l
II/-pll
pl(x)=
where the norm is
~
a~l"(x)
i=1
II/- p II= max(w Ix/I fix)
- p (x)1: x c x }
for a positive weight function w (x).
Characterization theorem [1] For each p (x) ~ K, define the following subsets of X : X+ 1 = (x E X:
S
can then be evaluated by interpolation for all
xEX. Execute the exchange algorithm with k = 1 and determine the subsets y2, U 2,L 2 as described below. Step k. Solve the n + 1 equations, i = 1. . . . . n + l :
W(X)[/(x)--p(x)] = II/-pll},
1-~-{x E X : w ( x ) [ / ( x ) - p ( x ) ]
=
-IIf-p[I},
X+2 = {x ~ X : p(x) = l(x)},
i6 yk,
i=1
~. aj~oj(xi k k ) = u (Xki), i ~ U k,
(1)
i=1
X _ 2 = {x E X : p ( x ) =/~ (x)},
Xp = X + I u X + 2 u X - 1 u X
k
O~i¢j(xikk ) _ _ ~ 1 ~/ O~n+X = f ( x k ) ,
2.
~.
k k ai~oi(xi)=l(x~),
i6L k
k=l
The theorem is then: p ~ K is the unique best constrained approximation to f if and only if there exist n + 1 consecutive points x~ < x 2 < • • •
for a ,k+ l . (The subset yk will be empty if n is too small, i.e. K = 0 ) . Define ek=[o~ k,+1[. The
F. Grenez / Constrained approximation for FIR filters
327
tions are
polynomial
p(x) = ~ aj cos 2 ( / ' - 1)~rx,
E
j=l
]=1
~j(x) = cos 2(/. - 1)wx, can be evaluated by interpolation. Execute the exchange algorithm and determine the subsets yk+l, uk+l, Lk+l (see below). The algorithm terminates when e k : IIf-p~[I, or, more practically, when the subsets yk, U k, L k remain unchanged by the exchange algorithm (on a dense grid replacing X).
and X is a compact subset of [0, 0.5]. Each extremal point x~ is characterized by its type tk: t/k = 0
if i~ yk,
tk = 1
if i c U k,
tk=-I
i f i ~ L k.
Exchange algorithm Each extremal point x k is replaced by the point k+l xi where the local maximum of yk occurs, with preservation of the alternation property of the characterization theorem, where: y/k = max{E/k, M/k, m/k},
Initially, at step 1, all extremal points are assigned to Y~: t~=O,
i=1 ..... n+l.
The system of equations (1) can be solved analytically f o r O~nk+l as: n+l
E k = max[w(x)]f(x)-pk(x)l-ek], Mki = max[pk (x ) - u(x )],
y~
k
:
(~ n + l
icYk+J
ify/k=E/k,
i 6 U k+l
ifyk=M/k,
i 6 L k+l
if y k = m k.
So, the new extremal point is the point where either the weighted error function is maximum, or the overstepping above u(x) is maximum, or the overstepping below l(x) is maximum.
3. Program implementation The classical program of Parks and McClellan for unconstrained Chebyshev approximation [3] can be converted into a general program for constrained approximation. The main modifications are described in this section. In the case of linear phase FIR filter design, the approximating func-
k
)
n+l
(2)
(--1)i-laki E i / W ( X k i )
m k = max[l(x ) -- p k(x)], and the type of the new extremal point is determined as follows:
i=1
k
aiq(xi
i=1
where n+l 1 a/k = l-I i=1 cos 2wx/k - c o s 2"trx~ '
i~i
~f(xki) if t/k=O, q(xk)=J'gU(Xi) if t / k = l , [ l ( x k) if t/k - 1 , {1 ek=_0
if t k =0, ift/k#0.
The denominator of (2) will vanish if all ek are zero, i.e. if the degree of the approximating polynomial is too small for the imposed constraints. The polynomial p k (X) can then be evaluated at any point x ~ X by the Lagrange interpolation formula:
pk(x ) _ i=1 i=1 Vol. 5, No. 4, July 1983
F. Grenez / Constrained approximation for FIR filters
328
where
Id.i =
A
cos 2"trx - cos 2rrx k'
A/k=fi 1 j=l cos 2wx/k - c o s 2wx k'
k
f(Xi)+(--l)~
cki =)u(xki) [I(x k)
k i O/n+1
if tk =0, if t k = 1, if t k = - 1 .
The exchange algorithm can be conducted in the same way as in the Parks and McClellan program. The local maximum of the function yk, with preservation of the sign of the error function associated with x k, is searched in a domain not exceeding the two adjacent extremal points Xk_l and xk+l. A further search is made for the endpoints of the domain. The algorithm terminates when the extremal points remain unchanged both in position and in type. The constrained Chebyshev approximation program is general for all types of FIR linear phase filters: case 2-3-4 filters (such as differentiators) are treated in the same way as in [3].
4. Applications
Linear-phase filters With the constrained Chebyshev approximation program, it is possible to define the allowable ripple in each band, and to obtain the best weighted Chebyshev approximation satisfying the constraints. The frequency response (approximating polynomial) will touch the imposed bounds u (x) and/or I(x) in some bands, and will be free in the other bands. A supplementary design parameter is provided by the weight function. Of course the weight function is only relevant where the Signal Processing
frequency response is not constrained by the bounds. The constrained Chebyshev approximation described in this paper can be compared to the constrained ripple design of [4] where the frequency response if forced to be tangent, at n + 1 points, to the upper and lower bounds whatever they could be. As opposed to the constrained Chebyshev approximation where the frequency response touches the bounds only when required by the optimality criterion.
Minimum-phase filters Although the constrained Chebyshev approximation program is concerned with linear-phase filters, it can be useful in the first step of the design of a minimum-phase filter when a linearphase prototype with double degree is needed. Indeed, the design of a minimum-phase filter is generally transposed to the approximation of the squared magnitude by a linear-phase prototype. The minimum-phase solution is then derived by a cepstral procedure [5, 6] or by spectral factorization [7]. As the squared magnitude must be nonnegative, the constrained Chebyshev approximation program can be applied with a lower bound l(x) = 0 in the stopbands and also in the transition bands if necessary (one-sided approximation). The double zeros of the prototype on the unit circle are known, as extremal points in the last iteration of the Remes algorithm. The methods proposed in [5, 7] to ensure the nonnegativity are only valid for low-pass or highpass filters. Indeed, in these methods, each band must necessarily be defined as either a passband (with free oscillation) or a stopband (with oscillation between zero and 8) and consequently the nonnegativity cannot be controlled in the transition bands. As the frequency response of a multiband filter (more than one transition band) is not guaranteed to be monotonic in the transition bands [8], this may lead to a negative squared magnitude in a transition band of a minimum phase filter.
< o_
BAND I 0000000000 .100000000 I000000000 1.000000000 I000000000 900000000 0~5110797 -3~.00~789974
C O M P U ] ER
] I I'-IF:
13450000E
~-01
EXIREMAL FR~OUENCIES .0148026( ) 0460526(+) 204934~(~) ~263158(+) 4000000(-) 4131579(+)
LOWER BAND EDgE UPPER BAND EDgE DESIRED VALUE WEIgHTINg UPPER B O U N D LOWER BOUND DEVIATION DEVIATION I N DB
3"/ = = = = = = = = = =
H( H( H( H( H( H( H( H( H( H(
.09~i053(+) .~7~3684(+) 5000000(+)
BAND 3 .200000000 .300000000 0.000000000 1.000000000 .010000000 0.000000000 .010000000 -40.000000000
37) 35) 33) 31) 29) 27) ~5) 23) 21) 19)
2) 4) 6) 8) i0) 1~) 14) l&) 18)
lb
.I000000( ) .~937500(-)
BAND 4 .300500000 .350000000 0.000000000 1.000000000 .030000000 0.000000000 .0~5110797 -3~.00~789974
H( H( H( H( H( H( H( H( H(
LgRID=
.1500000( .3~0~368(
) )
-.54463507E-03 -. 1767~724E-0~ -.~881597E-0~ -,~0~655~0E-0~ -.15b?~O~E-03 .93387948E-03 .10965919E-03 -.169&~74E-O~ -,29707806E-0~
BAND 5 .400000000 .500000000 1.000000000 1.000000000 1.020000000 .980000000 .0~0000000 -33.979400087
= = = = = = = = =
NgRID= ~45
Fig. 1. Output listing for a length 37 bandstop filter.
.0707237( ) 24934~1(-) 4707~37(-)
BAND ~ .150000000 .199500000 0.000000000 1.000000000 .030000000 0.000000000 .025110797 -32~00~789974
RESPONSE~*~*~ = .I077820&E-01 = .76301950E-03 = -.22783431E-01 = .~7545743E-02 = .45273228E-01 = -.2211019&E-0~ ~ -.93196749E-O1 = -.~1660005E-0~ = .3095~873E+00 = .50090647E+00
LPNO-[H =
~*IMPULSE H( I) H( 3) H( 5) H( 7) H( 9) H( I i ) H( ~ 3 ) H( ~ 5 ) H( 1 7 ) H( ] 9 )
FITTER
DANDPASS F I L T E R
F I N I T E IMPULSE RESPONSE LINEAR PHASE D I G I T A L DESIGN REMEZ EXCHANGE ALGORITHM WITH CONSTRAINTS
= = = = = = = = =
36) 34) 3~) 30) ~S) ~6) 24) ~) ~0)
.158~37(-) .3416184(-)
BAND
H( H( H( H( H( H( H( H( H(
~179b053( .3500000(
) )
cx
Ib
9
C~
F. Grenez / Constrained approximation for FIR filters
330
5. Examples The first example is a length 37 bandstop filter. The specifications are found on the output listing of Fig. 1 and the frequency response (or more precisely the zero-phase part of the frequency response) is given in Fig. 2. The desired value and
5ANOSFOP
FiLi'E~
riLTER LFNOTH:3/
..® .. I r--- . ~ /
1 ,,
'1
i ! ': ~
Ii
q- 13
o~
! I
/
I ®
®
:l' :
o ~, __.:
O~
L
,
1, lower constrained in bands 2-4, and totally constrained in bands 3-5. The extremal frequencies given in Fig. 1 are followed by (+) if the frequency response touches the upper bound (i U) or by (-) if the frequency response touches the lower bound (i ~ L). No weight function was used in this case (w(x) = 1 for all x ~X). The second example is a length 76 multiband filter (Fig. 3). A weight function, in inverse ratio to the allowable ripple, was used in this case. The frequency response is given in Fig. 4: only the upper bound is active in the passbands and the lower bound in the stopbands. Constraints have been added in the transition regions (bands 2 and 5) to suppress the parasitic ripples. The weight function was chosen small enough in bands 2-5 in order that the transition bands are only active by the upper and lower bounds, and not by the weighted error. Concerning the computation time, the constrained Chebyshev approximation program can be compared with the corresponding ParksMcClellan program with same filter length, grid density, desired value and weight function (but of course without upper and lower bounds). For the first and second example, the constrained approximation program requires respectively 5% and 9% more computer time than the Parks-McClellan program. In some cases, the increase in computer time can reach 50%, depending on the particular specifications.
!
.1
,2
S
,4
FRF@UFN~Y
Fig. 2. F r e q u e n c y r e s p o n s e for a l e n g t h 3 7 b a n d s t o p filter.
the lower bound are both zero in the stopbands whereas in the first passband the desired value and the upper bound are both equal to 1. The frequency response is upper constrained in band Signal Processing
6. Conclusions A general program for constrained Chebyshev approximation has been described. The program can be used for the design of any type of FIR linear phase filter, and also for the approximation of the squared magnitude of minimum phase filters.
F. Grenez / Constrained approximation ]:orFIR filters
+
* ,
~("J,d"~ ~i~O O~I;~,~" ~.~m~ ~ 0~01,~
*
+
~v~v OOOOOO~ ~OOOOOO~ OOOOOO~5~ OOOOOO0~ 0000000,81"ZO000C, lO~O" ~OOOOOOO~" m,-~OOOOOO~ '~'~OOOOO~-
• "od
'o
*
",n
*
I
I I I I I ,'~ ~ 0 0++t",l
~000000~'~-
I 0 ~
*
,0 r'l O N II e,,
*
IE
*
Z
* * *
. . . . .
~00]OOO~O~
~;~ ~ " ~ ,~"O ~ 01~ ,~;~ 15~ ~i'5 ~ O h. ,-~O ~ ~ ,-~~ ~ ~ ,.~i'5~ ~5~ ~0~ O~O~,'~ m~O~15~'~(5 ~i'~ ~ ' ~ ~ O ~ ' ~ O ~ ~ O --~In,-~OO OOhOh~,~h~h.,~.-~-~h~'lS~ •- ~ - - ~ 5 - - ~ - - ~ i " l ~ ' ; ~ O - - ~ i ~ O ~,-~(00.~,-~'~h.;~'0~(5~510~0 . . . . . . . . . . . . . . . . . . .
* =
~
OOOOOOOOOOOOOOOOOOO I I I I I I I i I I I I , I+l
I
OOOOOOOO ~OOOOOOOO OOOOOOOO OOOOOOOO ~OOOOOOOO Zl0~OOOOOO ~O~OOOOOO
mo,~oo * ,'~O~O~" 0~O~O~ * ,~0~010 * ~O~O~ * '~'~l,t,"~m * O~',t * *
~ 0 0 0 0 0 0
it II H II II II H II II II U II II II II II II II II ~ ~"~0m O ~ '~ ~ ~ O 0~"4"~ ~ O 0~" ~ m I- T T T T T 1 - ' r ' r - r ' r - r T T T T T T ' r
7^7^7 ,= *
000000,0;'-
~ v ~
OOOOF~O,I'~ 0000~0~ C,OOOO(50~0~ ZOOOO~OO~ '~OOOO,~O,-~,-~ ~OOOOO~ ~0~00000~5
• oo
o
*
* * *
~O~OOOO0J~¢
0 , 1 ~ ~1 ~ (~1(50.11~1 ,,.~OI ~ O.I ~ ~ ~ ,-~ 0 ,-~ 0
*
~-'~O~ ~)~O~ ~OO~
~O00000~l'~l
* * * * * ,
, * * *
....
OOOOOOm~ • OOOOO00]~ Io ~ , - ~ ;~~ ~ ,'~0~;~ ~ ~5,-~~ ~ ~ ~,-~0~ OOOOOO~
iIIIIIIIIIIIIIIIIII
331
~OIOO~ h.O,'~ O~ "001",O,0 0"01000" 0] ~ ~0.~~0 O,~ ~ 5 ~" . . . . .
"0 r,
* *
.g .=
* * *
=l
E
*
+~
*
000000~0~000000~ OOOOOO0]m OOOOOOh.~ c~OOOOOO,~ ~ZOOOOOO~ ~OOOOOO'~
+ + * ..... ~ 0~~ ~i'~ * 0~h-~ h-~0 * ~O~O~ * ~'O'~ * O ~"~0~"h- *
e-
*
m
,
m z ,..1
* ,
C
._=
, * ,
Z~ (.., "e
*
0]E
OOOOOOOO
I
* *
OOOOOOOO OOOOOOOO ~OOOOOOOO
m-~0J~O '~ 0~'~"m ~
* *
O O O O O O O O O O O O O O O O O O O ,~O0~OOOOOO I I I I I I I I I I I I I I I I I I+ m~O~O000 • W ~ L~ ~ b ~ l ~ l b ~ ~ l b ~ l ~ ~ 1 ~ L~I~ b ~ ~ OOOOOOOO •~ , ~ m ~ o o ~ m ~ o ~ m ~ . . . . . . . .
~ 0 ~ OOJ~m~ ~5~'~" O~i'~ ~"
* * * *
IIIIIIIIIIIIIIIIIII
6<0
,
~,,,
* *
01 f 3 .JlUZ
*
~
@ * * * * *
~ ~N ~ ZZ~ ~ ~
..... I I
R O 0 0 0 0 0 0 0
II II tl II II II II tl II II II It II It I1 II II It If
0 ~d ob
h-
~ ~ ~ 0 Z ~ ~
II • ~ ~
Z HI ~ ~
~
O
~
~
O
~
. . . . . . . . . . . . . . . . . . . II
I
II
?
III
Illllll[llllllllttllltlllllllltlllllll
~ ~
OOOOOO~
+
+
OOOOOO~ ~OOOOOO~ ZOOOO~O~ ~OOOOOOO~
vv~v~ ~ ~ O ~ ~
~ 0 ~ 0 0 0 0 0 ~
~
0000000~
E . . . . . . . . . . . . . vv . . . . ~ I I ~ I I I I ~ I I I I I I ! I I I I
,
*
oooooo~
~
0
I
mm UJWD Z Z > ~ Z Z ~ ZOO00 ~ 0 ~
0
~
~ W
*
* * * * * *
*
~d
o
L~ O O O O
m
Z
R ~ 0 ~ 0 ~
O
~
~ 0 ~ 0 ~ ~ 0 0 ~ 0 ~
*
* * *
Vol
5, No
4, July
1983
F. Orenez / Constrained approximation for FIR filters
332
References
~ULTIBAND FILTER FIL~E~ LENGTH=76
20
z
_111
c~
60 J J
'/J
0.
.l
2
.3
.4
so
.5
FREOUENC? Fig. 4. Frequency response for a length 76 multiband filter.
Signal Processing
[1] G.D. Taylor and M.J. Winter, "Calculation of best restricted ap[aroximations", SIAMJ. Numer. Anal., Vol. 7, No. 2, June 1970, pp. 248-255. [2] D.R. Gimlin, R.K. Cavin and M.C. Budge, "A multiple exchange algorithm for calculation of best restricted approximation", SIAM J. Numer. Anal., Vol. 11, No. 2, April 1974, pp. 219-231. [3] J.H. McClellan, T.W. Parks and L.R. Rabiner, "A computer program for designing optimum FIR linear phase digital filters, IEEE Trans., Vol. AU-21, No. 6, December 1973. [4] M.T. McCallig and B.J. Leon, "Constrained ripple design of FIR digital filters", IEEE Trans., Vol. CAS-25, No. 11, November 1978, pp. 893-902. [5] R. Boite and H. Leich, "A new procedure for the design of high order minimum phase FIR digital or CCD filters", SignalProcessing, Vol. 3, No. 3, 1981, pp. 101-108. [6] G.A. Mian and A.P, Nainer, "A fast procedure to design equiripple minimum-phase FIR Filters", IEEE Trans., Vol. CAS-29, No. 5, May 1982, pp. 327-331. [7] Y. Kamp and C.J. Wellekens, "Optimal design of minimum phase FIR filters", Proc. ICASSP-82, pp. 1817-1820. [8] L.R. Rabiner, J.F. Kaiser and R.W. Schafer, "Some considerations in the design of multiband finite-impulseresponse digital filters", IEEE Trans., Vol. ASSP-22, No. 6, December 1974, pp. 462-472.
325
D E S I G N OF L I N E A R OR M I N I M U M - P H A S E FIR FILTERS BY C O N S T R A I N E D C H E B Y S H E V A P P R O X I M A T I O N F. G R E N E Z , m e m b e r E U R A S I P
Service d'Electricit~ Ggndrale, C.P. 165, UniversitdLibre de Bruxelles, av. F. Roosevelt 50, 1050 Bruxelles, Belgium Received 15 December 1982 Revised 17 February 1983
Abstract. Constrained Chebyshev approximation is a known mathematical problem, and Remes' type algorithms for finding the solution have been described. It is shown in this paper that the classical Parks-McClellan program can be converted into a general program for constrained approximation. The program can serve for the design of both linear and minimum phase FIR digital filters.
Zusammeniassung.Die Tschebyshceff-Approximation mit Nebenbedigungen ist ein bekanntes mathematicsches Problem, welches mit Hilfe von Remes-artigen Algorithmen gelSst werden kann. Es wird in diesem Beitrag gezeigt, wie das klassische Program yon Parks und McClellan in einem allgemeinen Programm fiir die Approximation mit Nebenbedingungen modifiziert werden kann. Ein solches Programm kann zur Berechnung nichtrekursiver Filter mit linearer oder mit minimaler Phase gebraucht werden. R6sum6. L'approximation de Tchebycheff avec contraintes est un probl~me math6matique connu dont on peut obtenir la solution par des algorithmes de type Remez. On montre dans cet article comment le programme classique de Parks et McClellan peut ~tre converti en un programme g6n6ral d'approximation avec contraintes. Un tel programme peut servir pour le calcul des filtres non r6cursifs ~ phase lin6aire ou h phase minimale. Keywords.Constrained approximation, FIR filter.
1. Introduction
constraints of the form:
T h e C h e b y s h e v , or m i n i m a x , a p p r o x i m a t i o n of a c o n t i n u o u s f u n c t i o n f ( x ) by p o l y n o m i a l s is a classical p r o b l e m , a n d consists of finding the polyn o m i a l p ( x ) which m i n i m i z e s the u n i f o r m n o r m of the error f u n c t i o n :
[If-pll = max x
w(x)lf(x)-p(x)l.
T h e design of a linear phase F I R digital filter can be f o r m u l a t e d as a C h e b y s h e v a p p r o x i m a t i o n p r o b l e m a n d an efficient c o m p u t a t i o n p r o g r a m is available [3]. In some applications, it can be useful to restrict the class of a p p r o x i m a t i n g p o l y n o m i a l s by a d d i n g
l(x)<~p(x)<~u(x). F o r example, the allowable d e v i a t i o n in s o m e interval can be fixed to a given value, or the a p p r o x i m a t i n g p o l y n o m i a l can be r e q u i r e d to be nonnegative. This p r o b l e m is k n o w n in the literature as the c o n s t r a i n e d or restricted a p p r o x i m a t i o n . It is recalled here how the best restricted a p p r o x i m a t i o n can b e calculated a n d it is s h o w n that, for the design of F I R digital or C C D filters, the available p r o g r a m for u n c o n s t r a i n e d approxim a t i o n can b e modified to solve the m o r e g e n e r a l p r o b l e m of c o n s t r a i n e d a p p r o x i m a t i o n .
0165-1684/83/$3.00 © 1983, Elsevier Science Publishers B.V. (North-Holland)
F. Grenez / Constrained approximation for FIR filters
326
2. C o n s t r a i n e d C h e b y s h e v a p p r o x i m a t i o n
in Xp satisfying rr(xi) = (-1)'+Lo'(xl) where ~r(x) = - 1 if x ~ X 1 u X 2 and ~ r ( x ) = + l if x
Definitions
X + a k-) X + 2 .
Let X be a compact subset of [a, b ], qh • • • ~0, a Chebyshev set on X, and M the set of generalized approximating polynomials:
p (x) = ~ ~j~0;(x). i=l
Let f(x) be the function to be approximated and u (x), l(x) two auxiliary functions (called the upper and lower bound), with:
In other words, the best constrained approximation polynomial p (x) must exhibit n + 1 extremal points, where p(x) is alternatively above and below fix), and where either the weighted error function is m a x i m u m or where the polynomial touches one bound. This theorem generalizes the classical alternation theorem for unconstrained approximation.
Algorithm [2] l(x)
x6S. l(x)<~f(x)<~u(x), The subset K of M is defined as:
K = {p ~ M: l(x) <-p(x)<~ u (x), x ~ X}. K is the set of approximating polynomials whose ranges are restricted by I and u. Then the polynomial p*(x)~ K is the best constrained Chebyshev approximation to f if
II/-p*ll
= min p~K
A R e m e s ' type algorithm for the construction of the best constrained approximation can be described as follows. Step 1. Choose n + 1 extremai points x~ < x ~ < 1 • • • < x , + l and solve the n + 1 equations: 1
1
l
~ lai o~n+l - f ( x ] ) ,
) - , - , w(-U, )
/'=1
i=1 ..... n+l for an+l. 1 1 I. The polynomial Define e 1 = la.+l
II/-pll
pl(x)=
where the norm is
~
a~l"(x)
i=1
II/- p II= max(w Ix/I fix)
- p (x)1: x c x }
for a positive weight function w (x).
Characterization theorem [1] For each p (x) ~ K, define the following subsets of X : X+ 1 = (x E X:
S
can then be evaluated by interpolation for all
xEX. Execute the exchange algorithm with k = 1 and determine the subsets y2, U 2,L 2 as described below. Step k. Solve the n + 1 equations, i = 1. . . . . n + l :
W(X)[/(x)--p(x)] = II/-pll},
1-~-{x E X : w ( x ) [ / ( x ) - p ( x ) ]
=
-IIf-p[I},
X+2 = {x ~ X : p(x) = l(x)},
i6 yk,
i=1
~. aj~oj(xi k k ) = u (Xki), i ~ U k,
(1)
i=1
X _ 2 = {x E X : p ( x ) =/~ (x)},
Xp = X + I u X + 2 u X - 1 u X
k
O~i¢j(xikk ) _ _ ~ 1 ~/ O~n+X = f ( x k ) ,
2.
~.
k k ai~oi(xi)=l(x~),
i6L k
k=l
The theorem is then: p ~ K is the unique best constrained approximation to f if and only if there exist n + 1 consecutive points x~ < x 2 < • • •
for a ,k+ l . (The subset yk will be empty if n is too small, i.e. K = 0 ) . Define ek=[o~ k,+1[. The
F. Grenez / Constrained approximation for FIR filters
327
tions are
polynomial
p(x) = ~ aj cos 2 ( / ' - 1)~rx,
E
j=l
]=1
~j(x) = cos 2(/. - 1)wx, can be evaluated by interpolation. Execute the exchange algorithm and determine the subsets yk+l, uk+l, Lk+l (see below). The algorithm terminates when e k : IIf-p~[I, or, more practically, when the subsets yk, U k, L k remain unchanged by the exchange algorithm (on a dense grid replacing X).
and X is a compact subset of [0, 0.5]. Each extremal point x~ is characterized by its type tk: t/k = 0
if i~ yk,
tk = 1
if i c U k,
tk=-I
i f i ~ L k.
Exchange algorithm Each extremal point x k is replaced by the point k+l xi where the local maximum of yk occurs, with preservation of the alternation property of the characterization theorem, where: y/k = max{E/k, M/k, m/k},
Initially, at step 1, all extremal points are assigned to Y~: t~=O,
i=1 ..... n+l.
The system of equations (1) can be solved analytically f o r O~nk+l as: n+l
E k = max[w(x)]f(x)-pk(x)l-ek], Mki = max[pk (x ) - u(x )],
y~
k
:
(~ n + l
icYk+J
ify/k=E/k,
i 6 U k+l
ifyk=M/k,
i 6 L k+l
if y k = m k.
So, the new extremal point is the point where either the weighted error function is maximum, or the overstepping above u(x) is maximum, or the overstepping below l(x) is maximum.
3. Program implementation The classical program of Parks and McClellan for unconstrained Chebyshev approximation [3] can be converted into a general program for constrained approximation. The main modifications are described in this section. In the case of linear phase FIR filter design, the approximating func-
k
)
n+l
(2)
(--1)i-laki E i / W ( X k i )
m k = max[l(x ) -- p k(x)], and the type of the new extremal point is determined as follows:
i=1
k
aiq(xi
i=1
where n+l 1 a/k = l-I i=1 cos 2wx/k - c o s 2"trx~ '
i~i
~f(xki) if t/k=O, q(xk)=J'gU(Xi) if t / k = l , [ l ( x k) if t/k - 1 , {1 ek=_0
if t k =0, ift/k#0.
The denominator of (2) will vanish if all ek are zero, i.e. if the degree of the approximating polynomial is too small for the imposed constraints. The polynomial p k (X) can then be evaluated at any point x ~ X by the Lagrange interpolation formula:
pk(x ) _ i=1 i=1 Vol. 5, No. 4, July 1983
F. Grenez / Constrained approximation for FIR filters
328
where
Id.i =
A
cos 2"trx - cos 2rrx k'
A/k=fi 1 j=l cos 2wx/k - c o s 2wx k'
k
f(Xi)+(--l)~
cki =)u(xki) [I(x k)
k i O/n+1
if tk =0, if t k = 1, if t k = - 1 .
The exchange algorithm can be conducted in the same way as in the Parks and McClellan program. The local maximum of the function yk, with preservation of the sign of the error function associated with x k, is searched in a domain not exceeding the two adjacent extremal points Xk_l and xk+l. A further search is made for the endpoints of the domain. The algorithm terminates when the extremal points remain unchanged both in position and in type. The constrained Chebyshev approximation program is general for all types of FIR linear phase filters: case 2-3-4 filters (such as differentiators) are treated in the same way as in [3].
4. Applications
Linear-phase filters With the constrained Chebyshev approximation program, it is possible to define the allowable ripple in each band, and to obtain the best weighted Chebyshev approximation satisfying the constraints. The frequency response (approximating polynomial) will touch the imposed bounds u (x) and/or I(x) in some bands, and will be free in the other bands. A supplementary design parameter is provided by the weight function. Of course the weight function is only relevant where the Signal Processing
frequency response is not constrained by the bounds. The constrained Chebyshev approximation described in this paper can be compared to the constrained ripple design of [4] where the frequency response if forced to be tangent, at n + 1 points, to the upper and lower bounds whatever they could be. As opposed to the constrained Chebyshev approximation where the frequency response touches the bounds only when required by the optimality criterion.
Minimum-phase filters Although the constrained Chebyshev approximation program is concerned with linear-phase filters, it can be useful in the first step of the design of a minimum-phase filter when a linearphase prototype with double degree is needed. Indeed, the design of a minimum-phase filter is generally transposed to the approximation of the squared magnitude by a linear-phase prototype. The minimum-phase solution is then derived by a cepstral procedure [5, 6] or by spectral factorization [7]. As the squared magnitude must be nonnegative, the constrained Chebyshev approximation program can be applied with a lower bound l(x) = 0 in the stopbands and also in the transition bands if necessary (one-sided approximation). The double zeros of the prototype on the unit circle are known, as extremal points in the last iteration of the Remes algorithm. The methods proposed in [5, 7] to ensure the nonnegativity are only valid for low-pass or highpass filters. Indeed, in these methods, each band must necessarily be defined as either a passband (with free oscillation) or a stopband (with oscillation between zero and 8) and consequently the nonnegativity cannot be controlled in the transition bands. As the frequency response of a multiband filter (more than one transition band) is not guaranteed to be monotonic in the transition bands [8], this may lead to a negative squared magnitude in a transition band of a minimum phase filter.
< o_
BAND I 0000000000 .100000000 I000000000 1.000000000 I000000000 900000000 0~5110797 -3~.00~789974
C O M P U ] ER
] I I'-IF:
13450000E
~-01
EXIREMAL FR~OUENCIES .0148026( ) 0460526(+) 204934~(~) ~263158(+) 4000000(-) 4131579(+)
LOWER BAND EDgE UPPER BAND EDgE DESIRED VALUE WEIgHTINg UPPER B O U N D LOWER BOUND DEVIATION DEVIATION I N DB
3"/ = = = = = = = = = =
H( H( H( H( H( H( H( H( H( H(
.09~i053(+) .~7~3684(+) 5000000(+)
BAND 3 .200000000 .300000000 0.000000000 1.000000000 .010000000 0.000000000 .010000000 -40.000000000
37) 35) 33) 31) 29) 27) ~5) 23) 21) 19)
2) 4) 6) 8) i0) 1~) 14) l&) 18)
lb
.I000000( ) .~937500(-)
BAND 4 .300500000 .350000000 0.000000000 1.000000000 .030000000 0.000000000 .0~5110797 -3~.00~789974
H( H( H( H( H( H( H( H( H(
LgRID=
.1500000( .3~0~368(
) )
-.54463507E-03 -. 1767~724E-0~ -.~881597E-0~ -,~0~655~0E-0~ -.15b?~O~E-03 .93387948E-03 .10965919E-03 -.169&~74E-O~ -,29707806E-0~
BAND 5 .400000000 .500000000 1.000000000 1.000000000 1.020000000 .980000000 .0~0000000 -33.979400087
= = = = = = = = =
NgRID= ~45
Fig. 1. Output listing for a length 37 bandstop filter.
.0707237( ) 24934~1(-) 4707~37(-)
BAND ~ .150000000 .199500000 0.000000000 1.000000000 .030000000 0.000000000 .025110797 -32~00~789974
RESPONSE~*~*~ = .I077820&E-01 = .76301950E-03 = -.22783431E-01 = .~7545743E-02 = .45273228E-01 = -.2211019&E-0~ ~ -.93196749E-O1 = -.~1660005E-0~ = .3095~873E+00 = .50090647E+00
LPNO-[H =
~*IMPULSE H( I) H( 3) H( 5) H( 7) H( 9) H( I i ) H( ~ 3 ) H( ~ 5 ) H( 1 7 ) H( ] 9 )
FITTER
DANDPASS F I L T E R
F I N I T E IMPULSE RESPONSE LINEAR PHASE D I G I T A L DESIGN REMEZ EXCHANGE ALGORITHM WITH CONSTRAINTS
= = = = = = = = =
36) 34) 3~) 30) ~S) ~6) 24) ~) ~0)
.158~37(-) .3416184(-)
BAND
H( H( H( H( H( H( H( H( H(
~179b053( .3500000(
) )
cx
Ib
9
C~
F. Grenez / Constrained approximation for FIR filters
330
5. Examples The first example is a length 37 bandstop filter. The specifications are found on the output listing of Fig. 1 and the frequency response (or more precisely the zero-phase part of the frequency response) is given in Fig. 2. The desired value and
5ANOSFOP
FiLi'E~
riLTER LFNOTH:3/
..® .. I r--- . ~ /
1 ,,
'1
i ! ': ~
Ii
q- 13
o~
! I
/
I ®
®
:l' :
o ~, __.:
O~
L
,
1, lower constrained in bands 2-4, and totally constrained in bands 3-5. The extremal frequencies given in Fig. 1 are followed by (+) if the frequency response touches the upper bound (i U) or by (-) if the frequency response touches the lower bound (i ~ L). No weight function was used in this case (w(x) = 1 for all x ~X). The second example is a length 76 multiband filter (Fig. 3). A weight function, in inverse ratio to the allowable ripple, was used in this case. The frequency response is given in Fig. 4: only the upper bound is active in the passbands and the lower bound in the stopbands. Constraints have been added in the transition regions (bands 2 and 5) to suppress the parasitic ripples. The weight function was chosen small enough in bands 2-5 in order that the transition bands are only active by the upper and lower bounds, and not by the weighted error. Concerning the computation time, the constrained Chebyshev approximation program can be compared with the corresponding ParksMcClellan program with same filter length, grid density, desired value and weight function (but of course without upper and lower bounds). For the first and second example, the constrained approximation program requires respectively 5% and 9% more computer time than the Parks-McClellan program. In some cases, the increase in computer time can reach 50%, depending on the particular specifications.
!
.1
,2
S
,4
FRF@UFN~Y
Fig. 2. F r e q u e n c y r e s p o n s e for a l e n g t h 3 7 b a n d s t o p filter.
the lower bound are both zero in the stopbands whereas in the first passband the desired value and the upper bound are both equal to 1. The frequency response is upper constrained in band Signal Processing
6. Conclusions A general program for constrained Chebyshev approximation has been described. The program can be used for the design of any type of FIR linear phase filter, and also for the approximation of the squared magnitude of minimum phase filters.
F. Grenez / Constrained approximation ]:orFIR filters
+
* ,
~("J,d"~ ~i~O O~I;~,~" ~.~m~ ~ 0~01,~
*
+
~v~v OOOOOO~ ~OOOOOO~ OOOOOO~5~ OOOOOO0~ 0000000,81"ZO000C, lO~O" ~OOOOOOO~" m,-~OOOOOO~ '~'~OOOOO~-
• "od
'o
*
",n
*
I
I I I I I ,'~ ~ 0 0++t",l
~000000~'~-
I 0 ~
*
,0 r'l O N II e,,
*
IE
*
Z
* * *
. . . . .
~00]OOO~O~
~;~ ~ " ~ ,~"O ~ 01~ ,~;~ 15~ ~i'5 ~ O h. ,-~O ~ ~ ,-~~ ~ ~ ,.~i'5~ ~5~ ~0~ O~O~,'~ m~O~15~'~(5 ~i'~ ~ ' ~ ~ O ~ ' ~ O ~ ~ O --~In,-~OO OOhOh~,~h~h.,~.-~-~h~'lS~ •- ~ - - ~ 5 - - ~ - - ~ i " l ~ ' ; ~ O - - ~ i ~ O ~,-~(00.~,-~'~h.;~'0~(5~510~0 . . . . . . . . . . . . . . . . . . .
* =
~
OOOOOOOOOOOOOOOOOOO I I I I I I I i I I I I , I+l
I
OOOOOOOO ~OOOOOOOO OOOOOOOO OOOOOOOO ~OOOOOOOO Zl0~OOOOOO ~O~OOOOOO
mo,~oo * ,'~O~O~" 0~O~O~ * ,~0~010 * ~O~O~ * '~'~l,t,"~m * O~',t * *
~ 0 0 0 0 0 0
it II H II II II H II II II U II II II II II II II II ~ ~"~0m O ~ '~ ~ ~ O 0~"4"~ ~ O 0~" ~ m I- T T T T T 1 - ' r ' r - r ' r - r T T T T T T ' r
7^7^7 ,= *
000000,0;'-
~ v ~
OOOOF~O,I'~ 0000~0~ C,OOOO(50~0~ ZOOOO~OO~ '~OOOO,~O,-~,-~ ~OOOOO~ ~0~00000~5
• oo
o
*
* * *
~O~OOOO0J~¢
0 , 1 ~ ~1 ~ (~1(50.11~1 ,,.~OI ~ O.I ~ ~ ~ ,-~ 0 ,-~ 0
*
~-'~O~ ~)~O~ ~OO~
~O00000~l'~l
* * * * * ,
, * * *
....
OOOOOOm~ • OOOOO00]~ Io ~ , - ~ ;~~ ~ ,'~0~;~ ~ ~5,-~~ ~ ~ ~,-~0~ OOOOOO~
iIIIIIIIIIIIIIIIIII
331
~OIOO~ h.O,'~ O~ "001",O,0 0"01000" 0] ~ ~0.~~0 O,~ ~ 5 ~" . . . . .
"0 r,
* *
.g .=
* * *
=l
E
*
+~
*
000000~0~000000~ OOOOOO0]m OOOOOOh.~ c~OOOOOO,~ ~ZOOOOOO~ ~OOOOOO'~
+ + * ..... ~ 0~~ ~i'~ * 0~h-~ h-~0 * ~O~O~ * ~'O'~ * O ~"~0~"h- *
e-
*
m
,
m z ,..1
* ,
C
._=
, * ,
Z~ (.., "e
*
0]E
OOOOOOOO
I
* *
OOOOOOOO OOOOOOOO ~OOOOOOOO
m-~0J~O '~ 0~'~"m ~
* *
O O O O O O O O O O O O O O O O O O O ,~O0~OOOOOO I I I I I I I I I I I I I I I I I I+ m~O~O000 • W ~ L~ ~ b ~ l ~ l b ~ ~ l b ~ l ~ ~ 1 ~ L~I~ b ~ ~ OOOOOOOO •~ , ~ m ~ o o ~ m ~ o ~ m ~ . . . . . . . .
~ 0 ~ OOJ~m~ ~5~'~" O~i'~ ~"
* * * *
IIIIIIIIIIIIIIIIIII
6<0
,
~,,,
* *
01 f 3 .JlUZ
*
~
@ * * * * *
~ ~N ~ ZZ~ ~ ~
..... I I
R O 0 0 0 0 0 0 0
II II tl II II II II tl II II II It II It I1 II II It If
0 ~d ob
h-
~ ~ ~ 0 Z ~ ~
II • ~ ~
Z HI ~ ~
~
O
~
~
O
~
. . . . . . . . . . . . . . . . . . . II
I
II
?
III
Illllll[llllllllttllltlllllllltlllllll
~ ~
OOOOOO~
+
+
OOOOOO~ ~OOOOOO~ ZOOOO~O~ ~OOOOOOO~
vv~v~ ~ ~ O ~ ~
~ 0 ~ 0 0 0 0 0 ~
~
0000000~
E . . . . . . . . . . . . . vv . . . . ~ I I ~ I I I I ~ I I I I I I ! I I I I
,
*
oooooo~
~
0
I
mm UJWD Z Z > ~ Z Z ~ ZOO00 ~ 0 ~
0
~
~ W
*
* * * * * *
*
~d
o
L~ O O O O
m
Z
R ~ 0 ~ 0 ~
O
~
~ 0 ~ 0 ~ ~ 0 0 ~ 0 ~
*
* * *
Vol
5, No
4, July
1983
F. Orenez / Constrained approximation for FIR filters
332
References
~ULTIBAND FILTER FIL~E~ LENGTH=76
20
z
_111
c~
60 J J
'/J
0.
.l
2
.3
.4
so
.5
FREOUENC? Fig. 4. Frequency response for a length 76 multiband filter.
Signal Processing
[1] G.D. Taylor and M.J. Winter, "Calculation of best restricted ap[aroximations", SIAMJ. Numer. Anal., Vol. 7, No. 2, June 1970, pp. 248-255. [2] D.R. Gimlin, R.K. Cavin and M.C. Budge, "A multiple exchange algorithm for calculation of best restricted approximation", SIAM J. Numer. Anal., Vol. 11, No. 2, April 1974, pp. 219-231. [3] J.H. McClellan, T.W. Parks and L.R. Rabiner, "A computer program for designing optimum FIR linear phase digital filters, IEEE Trans., Vol. AU-21, No. 6, December 1973. [4] M.T. McCallig and B.J. Leon, "Constrained ripple design of FIR digital filters", IEEE Trans., Vol. CAS-25, No. 11, November 1978, pp. 893-902. [5] R. Boite and H. Leich, "A new procedure for the design of high order minimum phase FIR digital or CCD filters", SignalProcessing, Vol. 3, No. 3, 1981, pp. 101-108. [6] G.A. Mian and A.P, Nainer, "A fast procedure to design equiripple minimum-phase FIR Filters", IEEE Trans., Vol. CAS-29, No. 5, May 1982, pp. 327-331. [7] Y. Kamp and C.J. Wellekens, "Optimal design of minimum phase FIR filters", Proc. ICASSP-82, pp. 1817-1820. [8] L.R. Rabiner, J.F. Kaiser and R.W. Schafer, "Some considerations in the design of multiband finite-impulseresponse digital filters", IEEE Trans., Vol. ASSP-22, No. 6, December 1974, pp. 462-472.