EE32b Signals, Systems, and Transforms February 4, 2002
R. J. McEliece 162 Moore
Designing Discrete-Time Filters for Differentiation and Integration In these notes I will discuss the design and implementation of discrete-time filters for differentiation and integration. Together with the “sample, process in discrete time, and interpolate” paradigm of Chapter 7, these filters will allow us to build systems to differentiate and integrate continuous time signals using digital signal processing. 1. Discrete-Time Differentiation. As discussed in class (see also OW2, Sec. 7.4.1) the ideal digital differentiator has frequency response (1)
H(ejΩ ) = jΩ, for |Ω| ≤ π,
and repeats with period 2π. We wish to find a class of FIR filters which approximate this ideal as closely as possible. Since the impulse response corresponding to the ideal system (1) is real and odd,1 it is reasonable to investigate FIR approximations of the same form. Thus we will restrict our attention to filters characterized by system functions of the form (2)
N an n (z − z −n ), HN (z) = 2 n=1
and call N the order of the approximation. Substituting z = ejΩ in (2), we find that the frequency response of the N th order approximation is (3)
jΩ
HN (e ) = j
N
an sin(nΩ).
n=1
Comparing (1) to (3), we see that the design problem boils down to choosing the filter coefficients (an )N n=1 so that (4)
fN (Ω) =
N
an sin(nΩ) ≈ Ω,
n=1
In fact, this impulse response is h[0] = 0, h[n] = (−1)n /n, for n = 0. (See OW2, Ex. 7.2.) 1
1
in some appropriate sense. Although there are several possible criteria by which to interpret “≈,” for example minima; mean squared error, we will investigate filters which are optimal in the sense of maximal tangency to the ideal (1) at Ω = 0. If we use the Taylor series for sin x, viz., x5 x3 + − ··· 3! 5! x2m+1 (−1)m = (2m + 1)!
sin x = x −
m≥0
we can develop fN (Ω) as defined in (4) as a power series in Ω, viz., fN (Ω) = S1 Ω + S3 Ω3 + · · · , where the Sk ’s are linear combinations of the an ’s. Since fN (Ω) has N undetermined coefficients, it is reasonable to try to choose these coefficients so that S1 = 1, S3 = · · · = S2N −1 = 0, so that (5)
fN (Ω) = Ω + O(Ω2N +1 ).
For example, with N = 1, we have f1 (Ω) = a1 sin Ω = a1 (Ω −
Ω3 + · · ·), 3!
soo that if we choose a1 = 1, we get f1 (Ω) = Ω + O(Ω3 ), which is the closest match to the desired f (Ω) = Ω possible with just one degree of freedom. The system function for the corresponding filter is H1 (z) =
1 (z − z −1 ). 2
In the (discrete) time domain this filter is characterized by the difference equation y[n] =
1 (x[n + 1] − x[n − 1]). 2
which corresponds to the well-known (coninuous-time) approximation x (t) ≈
x(t + T ) − x(t − T ) . 2T 2
Similarly, with N = 2, we get Ω3 8Ω3 + · · ·) + a2 (2Ω − + · · ·) 3! 3! (a1 + 8a2 ) 3 = (a1 + 2a2 )Ω − Ω + ··· 3!
f1 (Ω) = a1 (Ω −
Thus if we select a1 andn a2 so that a1 + 2a2 = 1 a1 + 8a2 = 0, i.e., a1 = 4/3, a2 = −1/6, we get f2 (Ω) = Ω + O(Ω5 ). The corresponding system function is H2 (z) =
2 1 (z − z −1 ) − (z 2 − z −2 ), 3 12
with difference equation y[n] = −
1 2 2 1 x[n + 2] + x[n + 1] − x[n − 1] + x[n + 2]. 12 3 3 12
This corresponds to the not-so-well known approximation x (t) ≈
4 x(t + T ) − x(t − T ) 1 x(t + 2T ) − x(t − 2T ) − . 3 2T 3 4T
Figure 1 shows an implementation of a causal version of this system, viz., y[n] = −
1 2 2 1 x[n] + x[n − 1] − x[n − 3] + x[n − 4]. 12 3 3 12
y[n] -1/12
x[n]
2/3 z-1
1/12
-2/3 z-1
z-1
z-1
Figure 1. A causal FIR filter implementing a N = 2 differentiator.. 3
In general, we can expand fn (Ω) as a power series in Ω as follows: fN (Ω) =
N
an sin(nΩ)
n=1
=
N
an
n=1
(−1)m
m≥0
=
(−1)m
m≥0
(nΩ)2m+1 (2m + 1)!
N Ω2m+1 an n2m+1 (2m + 1)! n=1
Thus to ensure that (5) holds, we choose the filter coefficients (an ) so that N 1 if m = 0 2m+1 an n = 0 if m = 1 . . . , N − 1. n=1
This is a set of N linear equations in the N unknowns a1 , . . . , aN , which are easily solved by computer. Here is a short table of the solutions.
N N N N N
=1: =2: =3: =4: =5:
n=1
n=2
n=3
2/2 4/3 6/4 8/5 10/6
−1/6 −3/10 −6/15 −10/21
n=4
2/60 8/105 −1/140 20/168 −5/252
n=5
1/630
Figure 2 shows a plot of the optimiized functions f1 (Ω), . . . , f5 (Ω) vs. Ω. f (Ω) 3 N 2
N=5 1
N=1 -p
p - ÅÅÅÅ 2
p ÅÅÅÅ 2
-1 -2 -3
Figure 2. The first five “maximally flat” discrete-time differentiators.
4
p
Ω
2. Discrete-Time Integration. Motivated by the trapezoid rule and Simpson’s rule, perhaps, we consider discrete-time integrating filters characterized by difference equations of the form (6)
y[n] − y[n − N ] =
N
ak x[n − k],
k=0
where the filter coefficients (ak ) satisfy the symmetry condition (7)
ak = aN −k , for k = 0, . . . , N ,
and N is called the order of the filter. The system function corresponding to the difference equation (6) is N
(8) (9)
ak z −k 1 − z −N N ak z N/2−K = k=0 . z N/2 − z −N/2
HN (z) =
k=0
Our goal is to select the coefficients (ak ) so that the corresponding frequency response HN (ejΩ ) is as close as possible (in a suitable sense) to the ideal frequency response for a digital integrator, viz., (10)
H(ejΩ ) =
1 jΩ
as possible. For example, with N = 1, the system function in (9) is (recalling that a0 = a1 by (7)) H1 (z) =
a0 (z 1/2 + z −1/2 ) . z 1/2 − z −1/2
The corresponding frequency response is then easily seen to be 1 a0 cos(Ω/2) jΩ . (11) H1 (e ) = j sin(Ω/2) Comparing (10) to (11), we see that in this case the problem is to choose the coefficient a0 so that 1 a0 cos(Ω/2) ≈ , sin(Ω/2) Ω or alternatively (12)
Ωa0 cos(Ω/2) ≈ sin(Ω/2). 5
If we expand both sides of (12) as power series in Ω, using the Taylor series x5 x3 + − ··· 3! 5! x2 x4 cos x = 1 − + − ··· 2! 4! sin x = x −
(12) becomes (ignoring terms of order Ω5 and higher) a0 (Ω −
Ω3 Ω Ω3 + · · ·) ≈ − + ···. 8 2 48
Thus if we take as our measure of goodness the agreement of the power series at Ω = 0, the best choice for a0 is a0 = 1/2, in which case the filter’s system function is H1 (z) =
(1/2)(1 + z −1 ) , 1 − z −1
corresponding to the difference equation 1 y[n] = y[n − 1] + (x[n] + x[n − 1]), 2 which, as discussed in class, is the trapezoid rule for integration. Similarly, for N = 2 the system function is by (9) H2 (z) =
a1 + a0 (z + z −1 ) . z − z −1
Substituting z = ejΩ , this becomes 1 H2 (e ) = j jΩ
a1 /2 + a0 cos Ω sin Ω
.
Reasoning as in the N = 1 case, we see that the problem is to select the coefficients a0 and a1 so that (13)
Ω(a1 /2 + a0 cos Ω) ≈ sin Ω.
Expanding both sides of (13) in power series in Ω up to and including terms of order Ω3 , we get 1 (a1 /2 + a0 )Ω − (a0 /2)Ω3 + · · · ≈ Ω − Ω3 + · · · . 6 Thus we can get equality in (13) up to an including terms of order Ω4 if a0 and a1 are chosen so that 1 a0 + a1 = 1 2 1 1 = . a0 2 6 6
The solution is a0 = 1/3, a1 = 4/3, corresponding to the the system function H2 (z) =
1 3
+ 43 z −1 + 13 z −2 , 1 − z −2
with corresponding difference equation 1 4 1 y[n] = y[n − 2] + x[n] + x[n − 1] + x[n − 2], 3 3 3 which is, as discussed in class, nothing more or less than Simpson’s rule.
7
