Design of digital FIR notch filters

Design of digital FIR notch filters

S.C. Dutta Roy S.B. Jain B. Kumar

Indexing terms: Notch fillers, FIR digital filters

Abstract: A new procedure for designing digital FIR notch filters for a specified notch frequency wd and 3-dB rejection bandwidth Aco has been proposed. Design formulae for computing the required weights and the length N for the filter have been derived. Illustrative examples confirm- ing the new approach are also given.

1 Introduction

A notch filter highly attenuates a particular frequency component and leaves the others relatively unchanged.

Digital notch filters have many applications, such as the retrieval of sinusoids in noise, eliminating sinusoidal dis- turbances in biomedical, control and instrumentation systems, and tracking and enhancing time-varying narrowband/sinusoidal signals in wideband noise [1], among others. Different designs of IIR notch filters have been proposed by various researchers, including Carney [2] and Hirano [3]. In this paper we concern ourselves with the design of FIR notch filters. Tian Hu Yu and Mitra [4] used standard techniques for the design of FIR notch filters. Also, two techniques have been suggested whereby null width can be controlled L5J- In all these design methods, however, the weights or coefficients required for the filter structure are found by using computer-aided optimisation techniques. An analytical formulation of this problem, with exact mathematical for- mulae for the weights, does not appear to be available in the literature, and we here propose such an analytical design procedure for FIR notch filters. It has been shown that, through the proposed methodology, it is possible to realise exactly the desired notch frequency cod, besides meeting the specified rejection bandwidth Aco. We shall use the notation BW to indicate rejection bandwidth for notch filters in general.

Paper 1057K (E5, E10), first received 20th May 1993 and in revised form 5th January 1994

S.C. Dutta Roy is with the Department of Electrical Engineering, Indian Institute of Technology, Delhi, Hauz Khas, New Delhi 110016, India

S.B. Jain is with the Department of Electronics Engineering, Delhi College of Engineering, Kashmere Gate, Delhi 110006, India B. Kumar is with the Department of Electronics & Comm. Engng., Delhi Institute of Technology, Kashmere Gate, Delhi 110006, India

2 Design

Let | Hj(w) | be the magnitude response of an FIR digital notch filter of length N or order N — 1 (where N is assumed to be odd, in order to ensure integral delays in the realisation [6]), which meets the desired specifi- cations. We shall derive | Hd(w) | through the use of two maximally flat notch filters belonging to the class HJw), each of order N — 1 such that

(a) HJw) has m degrees of flatness at w = n (in the Butterworth sense), where m can assume (N — l)/2 differ- ent integral values;

(ft) the notch frequency of HJw) is wm, i.e. HJwm) = 0;and

(c) HJw) is positive for 0 ^ w < com and negative for com < a) ^ n.

Such an approach evades discontinuities of HJw) at co = wm and, as will be shown, it yields exact formulae for computation of the weights required for HJ^co). We shall choose m, the degree of flatness at co = n, in such a way that the resulting frequency response HJw) has a zero at to = vjm, where com is just short of cod. Let

HJw) = a s ico, n = •^{N - 1} (1)

where dt are the weights to be computed. Such an expres- sion for HJw) obviously represents the frequency response of an FIR linear-phase digital filter [7]. We impose the following optimality criteria for HJw):

flmM

d"HJw) dw"

d'HJto) dw"

, = 1 (2a)

: = - 1 (2ft)

= 0, u = 1, 2, ..., (2m - 1) (2c)

= 0, i; = 1, 2, ..., 2(n - m) + 1 (Id) It should be noted that the sum of the degrees of possible flatness of HJw) at co = 0 and w = n is N — 1. (For HJw) as chosen in eqn. 1, we require only n — 1 non- trivial equations in dt, besides the two constraints imposed by eqns. 2a and 2ft, to determine all the weights.

However, the odd-order derivatives (totalling n + 2) of HJw) are identically zero at co = 0 and n. This makes 2n{ = N — 1) derivatives of HJw), zero at w = 0 and it, i.e. HJw) has N — 1 degree of flatness reckoned at co = 0 and co = n together.) The integer m satisfies the relation 1 ^ m $ re. It is easily seen that, depending on the choice of n, the notch frequency wm can assume n discrete values

334 IEE Proc.-Vis. Image Signal Process., Vol. 141, No. 5, October 1994

a>1, 012, ..., a>n. Eqns. 2a-2d give the following n+ 1 nontrivial equations:

d0 + d, + d2 + • • • +

d0- dt + d2 +(-1)

0 + di - 22d2 + - - - + ( - l ) 0 - d , + 24d2 - • • • + ( - l )

0 +{~\)mdi + (-l)m+i2pd1 + ••• +(

0 + d, + 22d2 + • • • + 0 + d, + 24d2 + • • • +

d . = dB = n2dn = n*d. = nHn =

n

d] =

«4d = 1 - 1

0 0 0 0 0

0 2«d, = 0

where p = 2(m — 1); g = 2(n — m) and m = 2, 3, ..., n.

These equations can be put into matrix form as follows:

By using Crout's method [8], and following somewhat involved algebraic manipulations (see Appendix), eqn. 4 is transformed into the matrix equation

(5) where [aJJ is a triangular matrix with

0, i > j

J ~ ' ~ u

i = 0, 1, ..., m;

y = i + 1, i + 2, ..., n

*' //" + i~2l-4m + 1 V 2 m + / - 2

l A

b',=

i = m + 1, m + 2, ..., n;

j = i + 1, i + 2, . . . , n (6a) i = 0, 1, ..., m

- l ) i ' •») x 2 ( ' x

\ m - ly

i = m + 1, m + 2, ..., n (6b)

(3)

and

1, k = l 0, otherwise , = max. integer

s ! ( r - s ) !

(7a)

(7b)

(7c) Solving eqn. 5, the weighting coefficients d, are finally given by the following recursive formula:

i = n, n- 1, . . . , 0

(descending order) (8)

d, = b'

- £ a

j - i + I

where

i

£ ( • ) = () for / < fc

The values of the weights for H3(w), Hg(o>) and H15(w) when n = 15, for example, have been computed by using eqns. 6-8, and are given in Table 1. Proceeding in this way, we obtain the magnitude response curves \Hm(co)\

by varying m from 1 to n. It is observed that the 3 — dB BW of | Hm(w) | varies with m (keeping n constant). Fig. 1 shows the variation of BW with m for n = 9 and 15. It is Table 1: The weights d,, / = 0, 1, 2 15 torn = 15;m = 3. 8 and 15.

computed by using eqns. 6 and 8. The normalising factor is 22B

m = 3 m = 8

d0 0.1195582960D+09 rf, 0.2535146250D+09 d2 -0.1491262500D + 09 d3 0.4473787500D + 08 dA 0.1569750000D + 08 ds -0.2866363500D + 08 de 0.1726725000D+08 d-, -0.4148625000D + 07 da -0.1911000000D+07 da 0.2491125OOOD + 0 7 d, 0 -0.1365546000D + 07 rf,, 0.4845750000D + 06 d,2 -0.1183000000D + 06 d1 3 0.1942500000D + 05 dM -0.1950000000D + 04 d,. 0.9100000000D+02

0.O000O0O0O0D + 00 0.3312738000D + 09 0.0000000000D + 00 -0.85885800O0D + 08 0.0O0O0O0O00D + 00 0.3091888800D +08 0.O0O0O0O0O0D + 00 -0.1003860000D + 08 0.0O0O0O0O0OD + 00 0.2602600000D + 07 O.0O0O00O0O0D + 00 -0.4914000000D + 06 0.O0O0O0O0O0D + 00 0.5940000000D + 05 0.0000000000D + 00 -0.343200O0OOD + 04

-0.1908766960D + 09 0.1454226750D+09 0.1197598500D + 09 0.8649322500D + 08 0.5462730000D + 08 0.3004501500D + 08 0.1430715000D +08 0.5852925000D + 07 0.2035800000D + 07 0.5937750000D + 06 0.1425060000D + 06 0.2740500000D + 05 0.4060000000D + 04 0.4350000000D + 03 0.3000000000D + 02 0.1000000000D + 01 IEE Proc.-Vis. Image Signal Process., Vol. 141, No. 5, October 1994

also observed that the value of BW is maximum for m = l(n+ 1)/2J. Moreover, | Hm(<x>) | has the same BW for m = m0 and m = n + 1 — m0. This property of | Hm(o) \ facilitates the computation of weights corresponding to

0.64 0.54 0.53

0 51

10 12 U Fig. 1 Variation of 3-dB rejection bandwidth BW with m for n — 9 and 15

closest to, but less than, <od. The value of mi may be chosen from the n curves that Hm(w) assumes for different m. This, however, is a time-consuming procedure. We have therefore obtained the following empirical formula for m [:

m1 = integer part of [n(0.55 + 0.5 cos coj] (11) The empirical formulae given in eqns. 10 and 11 have been arrived at after modifying those proposed by Kaiser and Reed [9] and Hermann [10] in respect of maximally flat lowpass filters. These formulae have been verified for lengths up to 79.

A little reflection will show that, for Hm2(a>), where m2 = m, — 1, the notch frequency u>m2 is also closest to, but more thaji, cod. The desired notch frequency (a>d) is then obtained by linear combination of the frequency responses Hml(a>) and Hm2(o}). The values of weights D{ of the resulting filter are obtained as follows. Let

H((o) = a//ml(co) + pHm2(co) ' Dj cos ko where

com2 — o;m ] a»j — a>_,

(12a)

(12ft) (12c) m = n + 1 — m0 from those corresponding to m = m0 by

inspection using the relation

di mo = ( - l ) ' d , „ + ,_m o (9) Another significant observation about |ffm(a>)| is that BW progressively decreases as n increases. Fig. 2 shows the variation of BW for n varying from 3 to 37 and m = l(n + 1)/2J. This graph may now be used to find the value of optimum n required to obtain a specified BW for

Fig. 2 Variation of the rejection bandwidth BW with n for m = L(n + 1)12\ and n- 3,4, ...,37

| Hm(a>) |. An empirical formula which fits the curve shown in Fig. 2 is

n > integer H{(-R/B~W)2 - (n/BW) + 1}] (10) Thus, by using eqn. 10, we can readily compute the optimum n for a specified BW.

To arrive at the desired design for HJ,w), we first obtain Hml(w) such that its notch frequency {(oml) is

Dt = adj"1" + jidf2' (\2d) and

dlmj) = dt\m = mj (\2e) Eqns. 12ft and 12c follow from the facts that a and ji have to satisfy

a + ( 5 = 1 (13a) to ensure that H(0) = -H(n) = 1, and that the variations of Hml(w) and Hm2(w) near wd are approximately linear, so that

Hmi(w)

The resulting error due to the deviation from linearity of Hml(co) and Hm2(ct>) in the vicinity of at,, is very small and can be taken care of exactly, as discussed in the next Section.

The above-mentioned procedure yields a notch fre- quency of H(co) that is very close to the desired one (vjd).

However, it is observed that the 3-dB BW suffers from this procedure and becomes slightly higher than the specified one. This problem is overcome by increasing the value of n, found from eqn. 10, by 1, or equivalently, using the following revised empirical formula (instead of eqn. 10) for the computation of n:

n > integer l^{(n/BW)2 - (n/BW) + 3}] (14) The design procedure is formalised in the next Section.

3 Design procedure

To design an FIR notch filter HJ,CD) for the specified notch frequency (wd) and 3-dB rejection bandwidth (Aui), proceed as follows:

IEE Proc.-Vis. Image Signal Process., Vol. 141, No. 5, October 1994

(i) Choose the value_of_n from the empirical formula given by eqn. 14, using BW = Aco.

(ii) Find m^ from the formula given in eqn. 11, by using n found in step (i).

(iii) Use eqn. 8 to compute the weights d(/"j) corre- sponding to rtjj, = m, and m2 = ml — 1.

(iv) From the responses Hml(w) and Hm2(to), find the corresponding notch frequencies a>ml and o)m2, respec- tively. Compute a and p from eqns. 12b and 12c.

(v) Obtain weights D, by using eqn. I2d.

(vi) Find H(a>d) by using eqn. 12a. Ideally it should be zero, but because of deviations from the assumption leading to eqn. 12a, H((od) will be a small nonzero quan- tity £ of either sign. HJico) is then obtained as Hj(co) = H(o>) - c, which simply means that Do is modified to Do — s. This step may be called 'fine tuning', and helps to ensure that Hd(a>) has exactly zero value at w = <od. 4 Design examples

A computer program has been developed to evaluate the performance of the new design method. The design was carried out on a PC-AT/486 using double-precision arithmetic. Two example designs are presented here.

Example I: Let it be required to have a notch frequency (od = 2.5 radians with 3-dB rejection bandwidth Aco = 0.4 radians.

(i) The values of n and m1 are found to be 28 and 4, respectively, using eqns. 14 and 11.

(ii) The values of notch frequencies w^ and co3 are found to be 2.4064600 and 2.5164160 radians, respec- tively. Therefore, a = 0.149296 and /? = 0.8507039.

(iii) D, is obtained by linear combination of the weights corresponding to m, = 4 and m2 = 3 (i.e. <fj3)

and df). The value of the notch frequency for H((o) = YJi0 Di cos ico is 2.5007080 radians.

As H(w) does not have exactly the desired notch fre- quency, fine tuning is required. The value of H(2.5) is found to be 0.0023070 (£«,). We therefore modify Do to Do - Ei and obtain HJ&>) = H(a>) - £,. The notch fre- quency of HJtfo) is exactly 2.5 radians, as expected. The value of BW for HJco) is found to be 0.397 radians, which is better than the specified value of 0.4 radians.

The design therefore meets the given specifications.

5 Discussion

The design retains the maximal flatness (in the Butter- worth sense) of the passbands and achieves an exact null at the specified notch frequency. The exact mathematical

example no 2 example no '

Fig. 3 Example Example

1.0

Frequency responses \HmX(io)\, \Hml{a>)\ and \H/<o)\ =

!) + ^ Hw 2M — e \for examples 1 and 2

no. 1: «yd = 1 radian, 3-dB rejection bandwidth =0.375 radians no. 2: a>d = 2.5 radians, 3-dB rejection bandwidth = 0.4 radians

formulae for computing the weights needed constitute an attractive feature of the design. Also, this optimal design has been achieved without resorting to iterative tech- niques.

The rejection bandwidth {BW) can be made small if a sufficiently high value of n is chosen. For a given n, this filter provides a fixed range of notch frequencies varying from cod|? = n to cod\m = l. The range ((od\m=l - ci>d|m = .), however, increases with n, as shown in Fig. 4. The points

4 0

3 0 n 2 0

1 0 0

l i

• • | B

: i

\ \ . \

• \

; »

\ A

w

i NN^

n=39

n=9

i

, ' i

' 2 59j 1 1

B .f

i\ f '. :

'_ •

0 0 267 n/4

Example 2: Let it be desired to have a notch frequency a>d = 1 radian with 3-dB rejection bandwidth Aco = 0.375 radians.

The values of n and mt are found to be 32 and 26, respec- tively, using eqns. 14 and 11. The notch frequencies co26

and to25 are found to be 0.9427920 and 1.0172480 radians, respectively. The weights D{ corresponding to m, = 26 and m2 = 25 are calculated from eqn. 13 and the value of the notch frequency for H(w) = £ f jo D, cos ica is found to be 1.0002830 radians. For fine tuning, the value of f/(1.0) is found to be 0.0017113 ( = e2). We therefore modify Do to Do — E2 to arrive at the desired response Hj(co) for which the notch frequency is exactly equal to 1 radian. The value of BW for HJ^w) is found to be 0.373 radians, which is better than the specified value of 0.375 radians.

Fig. 3 shows the frequency responses \Hml(a>)\,

|ffm2(co)| and \HJ,w)\ for the two examples considered here. The frequency response H(CD) cannot be shown clearly on the graph, as it almost coincides with Hd(w).

IEE Proc.-Vis. Image Signal Process., Vol. 141, No. 5, October 1994

Fig. 4 The range of notch frequencies md for n varying from 3 to 39 on the left of co = n/2 on the graph shown in Fig. 4 corre- spond to the notch frequencies a>d\m=n, and those on the right are for to,,!„,«,. For example, with n = 9, wd\m = g = 0.55 radians and cod\m = i =2.59 radians, as shown by points A and A', respectively; and for n — 39, cod|m = 39 = 0.267 radians (point B) and coj|m = 1 = 2.88 radians (point B).

The following empirical inequality has been found to fit the range of notch frequencies for a given n:

21 — - + 0.51 U ( c o ,

cos"

2J--0.

A novel method has been suggested for designing an FIR notch filter. Formulae for computation of the required

weights have been derived and the procedure has been illustrated by design examples.

7 References

1 GLOVER, J.R.: 'Adaptive noise cancelling applied to sinusoidal interferences', IEEE Trans. Acoust. Speech, Signal Proc, 1977, 25, pp. 484-491

2 CARNEY, R.: 'Design of digital notch filter with tracking require- ments', IEEE Trans. Space Electron. Teiem., 1963, SET-9, pp. 109- 114

3 HIRANO, K., NISHIMURA, S., and MITRA, S.K.: 'Design of digital notch filters', IEEE Trans. Circuits Sys., 1974, CAS-21, pp.

540-546

4 YU, TIAN HU, MITRA, S.K., and BABIC HARVOJE: 'Design of linear phase notch filters', Sadhana, 1990, 15, (3), pp. 133-155 5 ER. M.H.: 'Designing notch filter with controlled null width', Signal

Proc, 1991, 24, pp. 319-329

6 RAB1NER, L.R., and GOLD, B.: 'Theory and application of digital signal processing' (Prentice Hall Inc., Englewood Cliffs, New Jersey, 1975)

7 OPPENHEIM, A.V., and SCHAFER, R.W.: 'Discrete-time signal processing' (Prentice Hall, Inc., Englewood Cliffs, NJ, 1989) 8 FROBERG, C: introduction to numerical analysis' (Reading, Mas-

sachusetts, Addison Wesley, 1969).

9 KAISER, J.F., and REED, W.A.:'Data smoothing using low pass digital filters', Rev. Sci. Instrum., 1977, 48, pp. 1447-1457 10 HERMANN, O.: 'On the approximation problem in nonrecursive

digital filter design', IEEE Trans. Circuit Theory, 1971, CT-18, pp.

411-413

11 RIORDON, J.: 'Combinatorial identities' (New York. John Wiley, 1968)

8 Appendix

Let [/I] = [ a , J , w h e r e i, j = 0, 1, 2, . . . , « ; [fc] = [fc0, bu

. . . , f tn]r; a n d T indicates t r a n s p o s e . I n C r o u t ' s m e t h o d , we t r a n s f o r m the elements of these m a t r i c e s to a'^ a n d b\, respectively, by using the following r e l a t i o n s [ 8 ] :

(16)

and b', = —

F o r e x a m p l e , t a k i n g n = 5 a n d m = 3, t h e m a t r i x e q u a - tion [af j][iij] = [ft,] is given by

1 1 0 0 0 0

1 - 1

1 - 1

1 1

1 1

— 22

2*

22

24

1 - 1

32

- 34

32

34

1 1 - 42

4*

42

4*

1 - 1

51

-54

51

5"

[do'

<*i

d2

d, dt

id,

1 - 1

0 0 0 0_

(18) Using eqns. 16 and 17, eqn. 18 is transformed to the fol- lowing:

1 1 0 0 0 _0

- -

- 1 2 1 1 1 1

1 0 - 4 16

4 16 1 1 1/22

1/2*

-3/26

3/27

I I - 2 - 4 8 16 112

As with Crout's method,

a'ij = a'ti =

X

Thus eqn.

1 1 1 D

' >j

' =i

1 0 4 - 4

64 640

we p u t

19 is modified to I

0 I

-

1 0 I 4 1 - 4 1

1 1 - 6

11

— 2 1

1 1 - 6

11

— 2 768_

<*i

d2

d3

d4

d.

do d1

d2 d}

di

~ 1 1 1/22

1/2*

- 3 / 26

3 / 27_ (19)

(20)

For general m and n, by using the combinatorial identity [10]

= (-1)1^{r + t - l (21)}

(17)

along with the method of induction, we finally arrive at eqns. 6a and 6b.

IEE Proc.-Vis. Image Signal Process., Vol. 141, No. S, October 1994