- Email: [email protected]
Direct design of recursive digital filters based on a new stability test
Direct Design of Recursive Digital Filters Based on a New Stability Test by V. RAMACHANDRAN
Department of Electrical Engineering, Concordia University, Montreal, Quebec, H3G 1M8, Canada C. S. GARGOUR
Ecole Polytechnique, University of Montreal, Montreal, Quebec, Canada M. AHMADI
and
M. T. H. BORAIE
Department of Electrical Engineering, University of Windsor, Wkdsor, Ontario, N9B 3P4, Canada
In this paper we present an alternative approach to the direct design of 1-D recursive digitalfilters satisfying prescribed magnitude specifications with or without constant group delay characteristic. This method uses an iterative method to calculate the coejkients of the$lter’s transfer,function and guarantees the stability of the designedjlter using a new stability test reported by Ramachandran and Gargour. To illustrate the usefulness of the technique, examples are given.
ABSTRACT:
I. Introduction There are two approaches to the design of recursive digital filters; the indirect and direct methods. In indirect methods, one of the classical analog-filter approximations (e.g. Butterworth, Bessel, Chebyshev, Elliptic) is used to generate an analog transfer function which is subsequently discretized by the use of the bilinear transformation method, the invariant-impulse response method (1). In direct techniques, on the other hand, the desired discrete-time transfer function is obtained directly from the given specifications, sometimes through the use of analytical methods (l), but more frequently through the use of an iterative method based on linear or nonlinear programming (2-S). The alternative method to be presented here is an iterative technique for the design of one-dimensional (1-D) recursive digital filter satisfying prespecified magnitude specification with or without constant group delay response. This method uses an unconstraint optimization method to approximate both the group delay and the magnitude response of the desired filter simultaneously if the constant group delay characteristic is required. The stability constraints which is based on a new stability test reported in (6) is imposed on the filter’s denominator polynomial through the variable substitution method described herein, hence ensuring the stability of the design filter.
ii-: The Frauklii Institute 0016+3032/84$3.00 +O.oO
407
V. Ramachandran
et al.
II. Characterization A one-dimensional function
of Recursive Filters recursive
digital
filter
D(z)
by its z-transfer
i$o mzi
N(z)
H,(z) = ~
is characterized
= ~
N izo 4i)z’
where N and D are polynomials in z = exp ST where s and T are complex variable and sampling period respectively. The design problem is to obtain the polynomial coefficients n(i) and d(i) such that the magnitude or the phase response of the filter, or both, approximate to the desired one beside maintaining the stability of the designed filter. The latter condition requires that D(z) # 0
IzI3 1.
(2)
The approximation can be carried out in z-domain (digital) by using any suitable nonlinear optimization method. The stability of the designed filter can be guaranteed by replacing all possible poles which are placed outside the unit circle by their mirror images with respect to the unit circle, hence, stabilizing the designed unstable filter. This stabilization procedure, however, has a drawback in which the phase response of the filter is not invariant under this stabilization process. In the method to be presented here this problem has been avoided through imposing necessary and sufficient conditions on the coefficients of a 1-D polynomial of any order having all its zeros inside the unit circle, based on a new stability test (6) throughout the optimization process.
ZZZ.Generation of 1-D Stable Polynomial It has been shown by Schussler coefficients as shown below
d(z) = dmzm+d,_Izm-l,..., = f.
(with d, positive) polynomial
D(z) can
of degree m having real
(7) that for a polynomial +dlz+d,
(3)
diz’
be decomposed
as the sum
of the mirror-image
F,(z) = +[D(z) + z”‘D(z - ‘)I and anti-mirror
(4)
image polynomial F,(z) = &D(z) - zmD(z- ‘)I.
For D(z) to have zeros inside the unit circle, the necessary are :
(5)
and sufficient conditions
(i) the zeros of F,(z) and F,(z) are located on the unit circle, Journalof 408
the Pergamon
Institute Press Ltd
Recursive
Digital
Filters
(ii) they are simple, (iii) they separate each other, and (iv)
II do
d,
< 1.
Ramachandran and Gargour (6) have used these results to generate a 1-D polynomial of order m which has all its zeros inside the unit circle using the following relationship : For m even D(z) = k, fi
(~~-2~cc~z+1)+(~~-1)“~~(z~--2~,z+1),
n=m 2
(6)
i=l
i=l
with k, > 1 and 1 > a, > p1 > ff2 > /I2 ... >/?,-I
> c(, > - 1.
(7)
For m odd D(z) = ko(z + 1) fi (z’ - 2aiz + 1) + (z - 1) fi i=l
(z’ - 28,~ + l),
i=l
n -
m-l 2
(8)
with k, > 1 and
(9)
l>a,>p,>a,>P,...>cc,>P,>-l and also in both cases ensure
f < 1 in Eq. (3). 1 m1
IV. Formulation of the Design Problem In the method to be presented here, Eq. (6) or Eq. (8) is assigned to the denominator of a 1-D transfer function Eq. (1). The error between the ideal and the designed magnitude response of the 1-D filter is calculated using the relationship E,(.kU
= IJ%Cexp
tiw,T)l I- IbCexp (_b,Vl I
(10)
where EM is the error of the magnitude response and lHrj and JH,I are the magnitude responses of the ideal and the designed filter respectively. In the same fashion the error between the ideal and the designed group delay response of the 1-D filter is calculated using the relationship &,,(j~,)
= r1 T- rwCexp
(k,T)l
(11)
where z1 is a constant representing the ideal group delay response of the filter and its value is chosen equal to the order of the filter (4) and r‘, is the group delay response of the designed filter. In this paper, two cases will be considered. (i) Formulation of the design problem for approximation of the magnitude response only. In this case we use the least mean square error criterion (I, norm) for Vol. 318, No. 6, pp. 407 413, December 1984 Printed in Great Britain
409
V. Ramachandran et al. the calculation
of the cost function
using Eq. (10) in the following
relationship
where I,, is the set of all discrete frequency points along the w axis in the passband and the stopband of the 1-D filter. Now any suitable nonlinear optimization technique can be used to calculate the parameters of the 1-D filter’s z-transfer function in such a way that El2 in Eq. (12) is minimized subject to the constraint of (7) or (9). This is a simple constraint optimization problem which can be solved either by utilization of any suitable optimization procedure with linear constraint, or to transform the problem to an unconstraint optimization problem by the use of the following variable substitution method (for m even) : ~1, = cos rc exp(-0:)
/?m_1 = cos
7~
exp[-(8:+@]
c&-r = cos rc exp [ - (0; + 13:+ 8,211
(13)
Similar conditions can be derived for m odd. Now any unconstraint optimization method can be used for the calculation of the new filter’s coefficients n and 8 to minimize E,, in Eq. (12). In this paper, we use the method of Fletcher and Powell (8) for this purpose. (ii) Formulation of the design problem for approximation of the magnitude and group delay response of the filter. In this case the general mean square error E, is calculated using Eqs. (lo), (11) in the following manner : (14) where I,, is the set of all discrete frequency point along the o axis in the passband and stopband and I, is the set of all discrete frequency points along the o axis in the passband of the 1-D filter respectively. Again, to design a 1-D filter satisfying prescribed magnitude and constant group delay response a, B, n should be calculated in such a way that E, in Eq. (14) is minimized subject to the constraint of (7) or (9). Using the steps exactly like those discussed in case (i) above, E, can be minimized either by utilization of any suitable nonlinear optimization routine with a linear constraint or to transform the problem to an unconstraint one using the variable substitution method previously discussed. V. Examples To illustrate the usefulness of the proposed method, several examples of recursive digital filters satisfying a prescribed magnitude response with or without constant group delay characteristics are given. In the first example, design of a bandpass filter with the following amplitude 410
Journal of the Frankhn Institute Pergamon Press Ltd.
Recursive
Digital
Filters
TABLEI
Values of the parameters of the designedjlter Denominator’s n, n, n2 n3 n4
specification
= 0.42153 = -0.01460 = -0.765959 = -0.016100 = 0.415030
xl =
coefficients
0.3867015
/I, = -0.0684533 t12 = -0.5061017
is required, f
0 for 0 d 101 < 0.5 rad/s
0 for 4.5 < 101 6 7
= 5 rad/s
where u, is considered to be 10 rad/s. The order of the filter is considered to be equal to four and a constant group delay characteristic is not required. Table I shows the values of n, a and j of the filter’s transfer function while Fig. 1 shows the magnitude response of the designed filter. In the second example, we design a lowpass filter with the following amplitude specification ( 1 for 0 < )o) $ 1
IH,(b4 = 011
-
IE
! 0 for i
2.5 d (01 < 7
= 5 rad/s
Olr
3
0607E
OO-
2
0506E
OO-
0405E
OO-
0304E
OO-
0203E
OO-
0.102E
OO-
0 609E-03 O.OOE 00
098E
00
0.20E
01
0.29E
01
0.39E
01
0.49E
01
w rad/s
FIG. 1. Amplitude
response
of the bandpass filter with the passband and w,,~ = 3 rad/s.
Vol. 318, No. 6, pp. 407 413, December 1984 Printed in Great Britain
edges of wPl = 2 rad/s
411
V. Ramachandran
et al. TABLE II Values of the parameters of the designedjilter Denominator’s n, n, n2 n3 n4
= 0.622769 = 1.031952 = 0.6625487 = -0.04637986 = -0.1994474
coefficients
cc1=
0.754979
/I1 =
0.3225849
tll = -0.8495334
and constant group delay characteristic. The order of the filter in this example is considered to be equal to four so we set z1 in Eq. (11) equal to four. Table II shows the value of the parameters of the designed filter while Figs. 2 and 3 show the magnitude and group delay response of the designed filter respectively.
VI. Conclusions In this paper we have presented an iterative method for the design of 1-D recursive digital filter satisfying a given magnitude response with or without constant group delay characteristics. Conditions are imposed on the parameters of the denominator polynomial of the filter’s z-transfer function to ensure the stability of the designed filter. These constraints are derived from a recently reported stability test given in (6). Then the constraint optimization problem of calculating the parameters of the designed filter
3418E
00
0 OOE
\
00
0.98E
00
X ml”
=
X,,,
-0.49000E
0.20E
O.OOOOOE
00 01
01
029E
01
039E
01
049EOl
Ymi, - 0.25792E -01 Ymai - 0 10054 E 01
FIG. 2. Amplitude response of the lowpass filter with the passband and stopband edges of
1 rad/s and 3 rad/s, respectively.
412
Journal
of the
Franklin Institute Pergamon Press Ltd.
Recursiw 0258E
Digital
Filters
02 r
0 236~ 0213E 0 19lE 0 169E 0 146E
0
IOZE
0 792E 0 569E 0 345E 0 122E 0 OOE 00
098E
00
OZOE w
FIG. 3. Group
is transformed
to an unconstraint
variable
substitution
designed
method.
method.
01
delay response
optimization Examples
Acknowledgement The authors are grateful for the support
0.29E
01
039E
01
0 49E
01
rod/s
of the lowpass filter.
problem
are given
by the application
to show
the usefulness
of the of the
of this project by the NSERC of Canada.
References (1) A. Antoniou, “Digital Filters: Analysis and Design”, McGraw Hill, New York, 1979. (2) A. T. Johnson, Jr., “Simultaneous magnitude and phase equalization using digital filters”, IEEE Trans., CAS-25, pp. 319-321, 1978. (3) G. C. Maenhout and W. Steenaart, “A direct approximation technique for digital filters and equalizers”, IEEE Trans., CT-20, pp. 548-555, 1973. (4) A. Chottera and G. A. Jullien, “A linear programming approach to recursive digital filter design with linear phase”, 1EEE Trans. Circuits Syst., Vol. CAS-29, No. 3, pp. 139-149, March 1982. (5) C. Charalambous and A. Antoniou, “Equalization ofrecursive digital filters”, IEEE Proc., Vol. 127, pt. G, No. 5, pp. 219-225, Oct. 1980. (6) V. Ramachandran and C. S. Gargour, “An implementation of a stability test of 1-D discrete system based on Schussler’s theorem and some consequent coefficient conditions”, J. Franklin Inst., Vol. 317, No. 5, pp. 341-358, May 1984. (7) H. W. Schussler, “A stability theorem for discrete systems”, IEEE Trans. Acoustics, Speech and Signal Processing, Vol. ASSP-24, No. 1, pp. 87789, Feb. 1976. (8) R. Fletcher and M. J. D. Powell, “A rapid descent method for minimization”, Computer J., Vol. 6, pp. 163-168, 1963.
Vol. 318, No 6, pp. 4.07-413, December Printed in Great Britam
IY84
413
Department of Electrical Engineering, Concordia University, Montreal, Quebec, H3G 1M8, Canada C. S. GARGOUR
Ecole Polytechnique, University of Montreal, Montreal, Quebec, Canada M. AHMADI
and
M. T. H. BORAIE
Department of Electrical Engineering, University of Windsor, Wkdsor, Ontario, N9B 3P4, Canada
In this paper we present an alternative approach to the direct design of 1-D recursive digitalfilters satisfying prescribed magnitude specifications with or without constant group delay characteristic. This method uses an iterative method to calculate the coejkients of the$lter’s transfer,function and guarantees the stability of the designedjlter using a new stability test reported by Ramachandran and Gargour. To illustrate the usefulness of the technique, examples are given.
ABSTRACT:
I. Introduction There are two approaches to the design of recursive digital filters; the indirect and direct methods. In indirect methods, one of the classical analog-filter approximations (e.g. Butterworth, Bessel, Chebyshev, Elliptic) is used to generate an analog transfer function which is subsequently discretized by the use of the bilinear transformation method, the invariant-impulse response method (1). In direct techniques, on the other hand, the desired discrete-time transfer function is obtained directly from the given specifications, sometimes through the use of analytical methods (l), but more frequently through the use of an iterative method based on linear or nonlinear programming (2-S). The alternative method to be presented here is an iterative technique for the design of one-dimensional (1-D) recursive digital filter satisfying prespecified magnitude specification with or without constant group delay response. This method uses an unconstraint optimization method to approximate both the group delay and the magnitude response of the desired filter simultaneously if the constant group delay characteristic is required. The stability constraints which is based on a new stability test reported in (6) is imposed on the filter’s denominator polynomial through the variable substitution method described herein, hence ensuring the stability of the design filter.
ii-: The Frauklii Institute 0016+3032/84$3.00 +O.oO
407
V. Ramachandran
et al.
II. Characterization A one-dimensional function
of Recursive Filters recursive
digital
filter
D(z)
by its z-transfer
i$o mzi
N(z)
H,(z) = ~
is characterized
= ~
N izo 4i)z’
where N and D are polynomials in z = exp ST where s and T are complex variable and sampling period respectively. The design problem is to obtain the polynomial coefficients n(i) and d(i) such that the magnitude or the phase response of the filter, or both, approximate to the desired one beside maintaining the stability of the designed filter. The latter condition requires that D(z) # 0
IzI3 1.
(2)
The approximation can be carried out in z-domain (digital) by using any suitable nonlinear optimization method. The stability of the designed filter can be guaranteed by replacing all possible poles which are placed outside the unit circle by their mirror images with respect to the unit circle, hence, stabilizing the designed unstable filter. This stabilization procedure, however, has a drawback in which the phase response of the filter is not invariant under this stabilization process. In the method to be presented here this problem has been avoided through imposing necessary and sufficient conditions on the coefficients of a 1-D polynomial of any order having all its zeros inside the unit circle, based on a new stability test (6) throughout the optimization process.
ZZZ.Generation of 1-D Stable Polynomial It has been shown by Schussler coefficients as shown below
d(z) = dmzm+d,_Izm-l,..., = f.
(with d, positive) polynomial
D(z) can
of degree m having real
(7) that for a polynomial +dlz+d,
(3)
diz’
be decomposed
as the sum
of the mirror-image
F,(z) = +[D(z) + z”‘D(z - ‘)I and anti-mirror
(4)
image polynomial F,(z) = &D(z) - zmD(z- ‘)I.
For D(z) to have zeros inside the unit circle, the necessary are :
(5)
and sufficient conditions
(i) the zeros of F,(z) and F,(z) are located on the unit circle, Journalof 408
the Pergamon
Institute Press Ltd
Recursive
Digital
Filters
(ii) they are simple, (iii) they separate each other, and (iv)
II do
d,
< 1.
Ramachandran and Gargour (6) have used these results to generate a 1-D polynomial of order m which has all its zeros inside the unit circle using the following relationship : For m even D(z) = k, fi
(~~-2~cc~z+1)+(~~-1)“~~(z~--2~,z+1),
n=m 2
(6)
i=l
i=l
with k, > 1 and 1 > a, > p1 > ff2 > /I2 ... >/?,-I
> c(, > - 1.
(7)
For m odd D(z) = ko(z + 1) fi (z’ - 2aiz + 1) + (z - 1) fi i=l
(z’ - 28,~ + l),
i=l
n -
m-l 2
(8)
with k, > 1 and
(9)
l>a,>p,>a,>P,...>cc,>P,>-l and also in both cases ensure
f < 1 in Eq. (3). 1 m1
IV. Formulation of the Design Problem In the method to be presented here, Eq. (6) or Eq. (8) is assigned to the denominator of a 1-D transfer function Eq. (1). The error between the ideal and the designed magnitude response of the 1-D filter is calculated using the relationship E,(.kU
= IJ%Cexp
tiw,T)l I- IbCexp (_b,Vl I
(10)
where EM is the error of the magnitude response and lHrj and JH,I are the magnitude responses of the ideal and the designed filter respectively. In the same fashion the error between the ideal and the designed group delay response of the 1-D filter is calculated using the relationship &,,(j~,)
= r1 T- rwCexp
(k,T)l
(11)
where z1 is a constant representing the ideal group delay response of the filter and its value is chosen equal to the order of the filter (4) and r‘, is the group delay response of the designed filter. In this paper, two cases will be considered. (i) Formulation of the design problem for approximation of the magnitude response only. In this case we use the least mean square error criterion (I, norm) for Vol. 318, No. 6, pp. 407 413, December 1984 Printed in Great Britain
409
V. Ramachandran et al. the calculation
of the cost function
using Eq. (10) in the following
relationship
where I,, is the set of all discrete frequency points along the w axis in the passband and the stopband of the 1-D filter. Now any suitable nonlinear optimization technique can be used to calculate the parameters of the 1-D filter’s z-transfer function in such a way that El2 in Eq. (12) is minimized subject to the constraint of (7) or (9). This is a simple constraint optimization problem which can be solved either by utilization of any suitable optimization procedure with linear constraint, or to transform the problem to an unconstraint optimization problem by the use of the following variable substitution method (for m even) : ~1, = cos rc exp(-0:)
/?m_1 = cos
7~
exp[-(8:+@]
c&-r = cos rc exp [ - (0; + 13:+ 8,211
(13)
Similar conditions can be derived for m odd. Now any unconstraint optimization method can be used for the calculation of the new filter’s coefficients n and 8 to minimize E,, in Eq. (12). In this paper, we use the method of Fletcher and Powell (8) for this purpose. (ii) Formulation of the design problem for approximation of the magnitude and group delay response of the filter. In this case the general mean square error E, is calculated using Eqs. (lo), (11) in the following manner : (14) where I,, is the set of all discrete frequency point along the o axis in the passband and stopband and I, is the set of all discrete frequency points along the o axis in the passband of the 1-D filter respectively. Again, to design a 1-D filter satisfying prescribed magnitude and constant group delay response a, B, n should be calculated in such a way that E, in Eq. (14) is minimized subject to the constraint of (7) or (9). Using the steps exactly like those discussed in case (i) above, E, can be minimized either by utilization of any suitable nonlinear optimization routine with a linear constraint or to transform the problem to an unconstraint one using the variable substitution method previously discussed. V. Examples To illustrate the usefulness of the proposed method, several examples of recursive digital filters satisfying a prescribed magnitude response with or without constant group delay characteristics are given. In the first example, design of a bandpass filter with the following amplitude 410
Journal of the Frankhn Institute Pergamon Press Ltd.
Recursive
Digital
Filters
TABLEI
Values of the parameters of the designedjlter Denominator’s n, n, n2 n3 n4
specification
= 0.42153 = -0.01460 = -0.765959 = -0.016100 = 0.415030
xl =
coefficients
0.3867015
/I, = -0.0684533 t12 = -0.5061017
is required, f
0 for 0 d 101 < 0.5 rad/s
0 for 4.5 < 101 6 7
= 5 rad/s
where u, is considered to be 10 rad/s. The order of the filter is considered to be equal to four and a constant group delay characteristic is not required. Table I shows the values of n, a and j of the filter’s transfer function while Fig. 1 shows the magnitude response of the designed filter. In the second example, we design a lowpass filter with the following amplitude specification ( 1 for 0 < )o) $ 1
IH,(b4 = 011
-
IE
! 0 for i
2.5 d (01 < 7
= 5 rad/s
Olr
3
0607E
OO-
2
0506E
OO-
0405E
OO-
0304E
OO-
0203E
OO-
0.102E
OO-
0 609E-03 O.OOE 00
098E
00
0.20E
01
0.29E
01
0.39E
01
0.49E
01
w rad/s
FIG. 1. Amplitude
response
of the bandpass filter with the passband and w,,~ = 3 rad/s.
Vol. 318, No. 6, pp. 407 413, December 1984 Printed in Great Britain
edges of wPl = 2 rad/s
411
V. Ramachandran
et al. TABLE II Values of the parameters of the designedjilter Denominator’s n, n, n2 n3 n4
= 0.622769 = 1.031952 = 0.6625487 = -0.04637986 = -0.1994474
coefficients
cc1=
0.754979
/I1 =
0.3225849
tll = -0.8495334
and constant group delay characteristic. The order of the filter in this example is considered to be equal to four so we set z1 in Eq. (11) equal to four. Table II shows the value of the parameters of the designed filter while Figs. 2 and 3 show the magnitude and group delay response of the designed filter respectively.
VI. Conclusions In this paper we have presented an iterative method for the design of 1-D recursive digital filter satisfying a given magnitude response with or without constant group delay characteristics. Conditions are imposed on the parameters of the denominator polynomial of the filter’s z-transfer function to ensure the stability of the designed filter. These constraints are derived from a recently reported stability test given in (6). Then the constraint optimization problem of calculating the parameters of the designed filter
3418E
00
0 OOE
\
00
0.98E
00
X ml”
=
X,,,
-0.49000E
0.20E
O.OOOOOE
00 01
01
029E
01
039E
01
049EOl
Ymi, - 0.25792E -01 Ymai - 0 10054 E 01
FIG. 2. Amplitude response of the lowpass filter with the passband and stopband edges of
1 rad/s and 3 rad/s, respectively.
412
Journal
of the
Franklin Institute Pergamon Press Ltd.
Recursiw 0258E
Digital
Filters
02 r
0 236~ 0213E 0 19lE 0 169E 0 146E
0
IOZE
0 792E 0 569E 0 345E 0 122E 0 OOE 00
098E
00
OZOE w
FIG. 3. Group
is transformed
to an unconstraint
variable
substitution
designed
method.
method.
01
delay response
optimization Examples
Acknowledgement The authors are grateful for the support
0.29E
01
039E
01
0 49E
01
rod/s
of the lowpass filter.
problem
are given
by the application
to show
the usefulness
of the of the
of this project by the NSERC of Canada.
References (1) A. Antoniou, “Digital Filters: Analysis and Design”, McGraw Hill, New York, 1979. (2) A. T. Johnson, Jr., “Simultaneous magnitude and phase equalization using digital filters”, IEEE Trans., CAS-25, pp. 319-321, 1978. (3) G. C. Maenhout and W. Steenaart, “A direct approximation technique for digital filters and equalizers”, IEEE Trans., CT-20, pp. 548-555, 1973. (4) A. Chottera and G. A. Jullien, “A linear programming approach to recursive digital filter design with linear phase”, 1EEE Trans. Circuits Syst., Vol. CAS-29, No. 3, pp. 139-149, March 1982. (5) C. Charalambous and A. Antoniou, “Equalization ofrecursive digital filters”, IEEE Proc., Vol. 127, pt. G, No. 5, pp. 219-225, Oct. 1980. (6) V. Ramachandran and C. S. Gargour, “An implementation of a stability test of 1-D discrete system based on Schussler’s theorem and some consequent coefficient conditions”, J. Franklin Inst., Vol. 317, No. 5, pp. 341-358, May 1984. (7) H. W. Schussler, “A stability theorem for discrete systems”, IEEE Trans. Acoustics, Speech and Signal Processing, Vol. ASSP-24, No. 1, pp. 87789, Feb. 1976. (8) R. Fletcher and M. J. D. Powell, “A rapid descent method for minimization”, Computer J., Vol. 6, pp. 163-168, 1963.
Vol. 318, No 6, pp. 4.07-413, December Printed in Great Britam
IY84
413