1
2
3
frequency (radians/sample)
Figure 9.11: Magnitude frequency response in decibels of a 64-point lowpass FIR
filter designed with the Parks-McClellan algorithm.
Lee & Varaiya, Signals and Systems
393
9.4. FINITE IMPULSE RESPONSE (FIR) FILTERS
Probing Further: Decibels
The term “decibel” is literally one tenth of a bel, which is named after Alexander Gra-
ham Bell. This unit of measure was originally developed by telephone engineers at Bell
Telephone Labs to designate the ratio of the power of two signals.
Power is a measure of energy dissipation (work done) per unit time. It is measured in
watts for electronic systems. One bel is defined to be a factor of 10 in power. Thus, a
1000 watt hair dryer dissipates 1 bel, or 10 dB, more power than a 100 watt light bulb.
Let p1 = 1000 watts be the power of the hair dryer and p2 = 100 be the power of the
light bulb. Then the ratio is
log (
10 p1/ p2) = 1 bel,
or
10 log (
10 p1/ p2) = 10 dB.
Comparing against (9.24) we notice a discrepancy. There, the multiplying factor is
20, not 10. That is because the ratio in (9.24) is a ratio of amplitudes, not powers.
In electronic circuits, if an amplitude represents the voltage across a resistor, then the
power dissipated by the resistor is proportional to the square of the amplitude. Let a1
and a2 be two such amplitudes. Then the ratio of their powers is
10 log (
/ ) =
(
10 a2
1 a2
2
20 log10 a1/a2).
Hence the multiplying factor of 20 instead of 10. A 3 dB power ratio amounts to a
√
factor of 2 in power. In amplitudes, this is a ratio of
2. The edge of the passband of
a bandpass filter is often defined to be the frequency at which the power drops to half,
which is the frequency where the gain is 3 dB below the passband gain. The magnitude
√
frequency response here is 1/ 2 times the passband gain.
In audio, decibels are used to measure sound pressure. A statement like “a jet engine
at 10 meters produces 120 dB of sound,” by convention, compares sound pressure to
a defined reference of 20 micropascals, where a pascal is a pressure of 1 newton per
square meter. For most people, this is approximately the threshold of hearing at 1 kHz.
Thus, a sound at 0 dB is barely audible. A sound at 10 dB has 10 times the power. A
sound at 100 dB has 1010 times the power. The jet engine, therefore, would probably
make you deaf without ear protection.
394
Lee & Varaiya, Signals and Systems
9. FILTERING
9.5
Infinite impulse response (IIR) filters
The equation for an FIR filter (9.22) is a difference equation relating the output at index
n to the inputs at indices n − L + 1 through n. A more general form of this difference
equation includes more indices of the output,
L−1
M−1
∀ n ∈ Z, y(n) = ∑ x(n−m)b(m)+ ∑ y(n−m)a(m)
(9.25)
m=0
m=1
where L and M are positive integers, and b(m) and a(m) are filter coefficients, which are
usually real valued. If all the a coefficients are zero, then this reduces to an FIR filter
with impulse response b = h. However, if any a coefficient is non-zero, then the impulse
response of this filter will never completely die out, since the output keeps getting recycled
to affect future outputs. Such a filter is called an infinite impulse response or IIR filter, or sometimes a recursive filter.
Example 9.16: Consider a causal LTI system defined by the difference equation
∀ n ∈ Z, y(n) = x(n) + 0.9y(n − 1).
We can find the impulse response by letting x = δ, the Kronecker delta function.
Since the system is causal, we know that h(n) = 0 for n < 0 (see page 382). Thus,
y(n)
=
0
if n < 0
y(0)
=
1
y(1)
=
0.9
y(2)
=
0.92
y(3)
=
0.93
Noticing the pattern, we conclude that
y(n) = (0.9)nu(n),
where u is the unit step, seen before in (2.16),
1
if n ≥ 0
u(n) =
(9.26)
0
otherwise
The output y and the magnitude of the frequency response (in dB) are plotted in
Figure 9.12. Notice that this filter has about 20 dB of gain at DC, dropping to about
-6 dB at π radians/sample.
Lee & Varaiya, Signals and Systems
395
9.5. INFINITE IMPULSE RESPONSE (IIR) FILTERS
9.5.1
Designing IIR filters
Designing IIR filters amounts to choosing the a and b coefficients for (9.25). As with
FIR filters, how to choose these coefficients is a well-studied subject, and the results of
this study are (for the most part) available in easy-to-use software. There are four widely
used methods for calculating these coefficients, resulting in four types of filters called
Butterworth, Chebyshev 1, Chebyshev 2, and elliptic filters.
Example 9.17:
We can design one of each of the four types of filters using the
Matlab commands butter, cheby1, cheby2, ellip. The arguments for each
of these specify either a cutoff frequency, which is the frequency at which the
filter achieves -3 dB gain, or in the case of cheby2, the edge of the stopband. For
example, to get lowpass filters with a gain of about 1 from 0 to π/8 radians/sample,
and a stopband at higher frequencies, we can use the following commands:
N = 5;
Wn = 0.125;
[B1, A1] = butter(N, Wn);
[B2, A2] = cheby1(N,1,Wn);
[B3, A3] = cheby2(N,70,0.25);
[B4, A4] = ellip(N,1,70,Wn);
The returned values are vectors containing the a and b coefficients of (9.25). The
magnitude frequency response of the resulting filters is shown in Figure 9.13. In
that figure, we show only positive frequencies, since the magnitude frequency re-
sponse is symmetric. We also only show the frequency range from 0 to π, since the
frequency response is periodic with period 2π. Notice that the Butterworth filter has
the most gradual rolloff from the passband to the stopband. Thus, for a given filter
order, a Butterworth filter yields the weakest lowpass characteristic of the four. The
elliptic filter has the sharpest rolloff. The Butterworth filter, on the other hand, has
the smoothest characteristic. The Chebyshev 1 filter has ripple in the passband,
the Chebyshev 2 filter has ripple in the stopband, and the elliptic filter has ripple in
both.
In the above Matlab commands, N is the filter order, equal to L and M in (9.25).
It is a constraint of these filter design techniques that L = M in (9.25). Wn is the
cutoff frequency, as a fraction of π radians/sample. A cutoff frequency of 0.125
396
Lee & Varaiya, Signals and Systems
9. FILTERING
1
0.8
0.6
0.4
0.2
0
−5
0
5
10
15
20
25
30
discrete time index
30
20
10
amplitude (in dB)
0
−10 −3
−2
−1
0
1
2
3
frequency (radians/sample)
Figure 9.12: Impulse response (top) and magnitude frequency response in dB
(bottom) of a simple causal IIR filter.
Lee & Varaiya, Signals and Systems
397
9.6. IMPLEMENTATION OF FILTERS
means 0.125π = π/8 radians/sample. The “1” argument in the cheby1 command
specifies the amount of passband ripple that we are willing to tolerate (in dB).
The 70 in the cheby2 and ellip commands specifies the amount of stopband
attenuation that we want (in dB). Finally, the 0.25 in the cheby2 line specifies the
edge of the stopband.
9.6
Implementation of filters
We now have several ways to describe an LTI system (a filter). We can give a state-space
description, a frequency response, an impulse response, or a difference equation such as
(9.25) or (9.22). All are useful. The state-space description and the difference equations prove the most useful when constructing the filter.
A realization of a filter in hardware or software is called an implementation. Do not
confuse filter design with filter implementation. The term “filter design” is used in the
community to refer to the choice of frequency response, impulse response, or coefficients
for a difference equation, not to how the frequency response is implemented in hardware
or software. In this section, we talk about implementation.
The output y of an FIR filter with impulse response h and input x is given by (9.22).
The output of an IIR filter with coefficients a and b is given by (9.25). Each of these
difference equations defines a procedure for calculating the output given the input. We
discuss various ways of implementing these procedures.
9.6.1
Matlab implementation
If x is finite, and we can interpret it as an infinite signal that is zero outside the specified
range. Then we can compute the output of an FIR filter using Matlab’s conv function,
and the output of an IIR filter using filter.
Example 9.18: Consider an FIR filter with an impulse response of length L. If x
is a vector containing the P values of the input, and h is a vector containing the L
values of the impulse response, then
398
Lee & Varaiya, Signals and Systems
9. FILTERING
20
20
0
0
−20
−20
−40
−40
−60
amplitude (dB)
−60
amplitude (dB)
−80
−80
−100
−100
0
1
2
3
0
1
2
3
frequency (radians/sample)
frequency (radians/sample)
20
20
0
0
−20
−20
−40
−40
−60
−60
amplitude (dB)
amplitude (dB)
−80
−80
−100
−100
0
1
2
3
0
1
2
3
frequency (radians/sample)
frequency (radians/sample)
Figure 9.13: Four frequency responses for 5-th order IIR filters of type Butter-
worth (upper left), Chebyshev 1 (upper right), Chebyshev 2 (lower left), and Ellip-
tic (lower right).
Lee & Varaiya, Signals and Systems
399
9.6. IMPLEMENTATION OF FILTERS
y = conv(x, h);
yields a vector containing L + P − 1 values of the output.
This strategy, of course, only works for finite input data, since Matlab requires that the
input be available in a finite vector.
Discrete-time filters can be implemented using standard programming languages and us-
ing assembly languages on specialized processors (see boxes). These implementations do
not have the limitation that the input be finite.
9.6.2
Signal flow graphs
We can describe the computations in a discrete-time filter using a block diagram with
three visual elements, a unit delay, a multiplier, and an adder. In the convolution sum for
an FIR filter,
L−1
y(n) = ∑ h(m)x(n − m),
m=0
notice that at each n we need access to x(n), x(n − 1), x(n − 2), · · · , x(n − L + 1). We can
maintain this set of values by cascading a set of unit delay elements to form a delay line,
as shown in figure 9.14.
For each integer n, the output sample is the values in the delay line scaled by h(0), h(1),
···, h(L−1). To obtain the values in the delay line, we simply tap the delay line, as shown
in Figure 9.15. The triangular boxes denote multipliers that multiply by a constant (h(m),
in this case). The circles with the plus signs denote adders. The structure in Figure 9.15
is called a tapped delay line description of an FIR filter.
Diagrams of the type shown in Figure 9.15 are called signal flow graphs because they
describe computations by focusing on the flow of signals between operators. Signal flow
graphs can be quite literal descriptions of digital hardware that implements a filter. But
what they really describe is the computation, irrespective of whether the implementation
is in hardware or software.
An IIR filter can also be described using a signal flow graph. Consider the difference
equation in (9.25). A signal flow graph description of this equation is shown in Figure
9.16. Notice that the left side is structurally equivalent to the tapped delay line, but stood
400
Lee & Varaiya, Signals and Systems
9. FILTERING
Probing Further: Java implementation of an FIR filter
The following Java class implements an FIR filter:
1
class FIR {
2
private int length, count = 0;
3
private double[] delayLine, impResp;
4
FIR(double[] coefs) {
5
length = coefs.length;
6
impResp = coefs;
7
delayLine = new double[length];
8
}
9
double getOutputSample(double inputSample) {
10
delayLine[count] = inputSample;
11
double result = 0.0;
12
int index = count;
13
for (int i=0; i<length; i++) {
14
result += impResp[i] * delayLine[index--];
15
if (index < 0) index = length-1;
16
}
17
if (++count >= length) count = 0;
18
return result;
19
}
20
}
A class in Java has both data members and methods. The methods are procedures
that operate on the data, and may or may not take arguments or return values. In this
case, there are two procedures, “FIR” and “getOutputSample.” The first, lines 4-8,
is a constructor, which is a procedure that is called to create an FIR filter. It takes
one argument, an array of double-precision floating-point numbers that specify the
impulse response of the filter. The second, lines 9-19, takes a new input sample as an
argument and returns a new output sample. It also updates the delay line using circular
buffering, where the “count” member is used to keep track of where each new input
sample should go. It gets incremented (line 17) each time the getOutputSample()
method is called. When it exceeds the length of the buffer, it gets reset to zero. At all
times, it contains the L most recently received input data samples. The most recent one
is at index count in the buffer. The second most recent is at count - 1, or if that
is negative, at length - 1. Line 15 makes sure that the variable index remains
within the buffer as we iterate through the loop.
Lee & Varaiya, Signals and Systems
401
9.6. IMPLEMENTATION OF FILTERS
Probing Further: FIR filter in a programmable DSP
The following section of code is the assembly language for a programmable DSP, which
is a specialized microprocessor designed to implement signal processing functions ef-
ficiently in embedded systems (such as cellular telephones, digital cordless telephones,
digital audio systems, etc.). This particular code is for the Motorola DSP56000 family
of processors.
1
fir
movep
x:input,x:(r0)
2
clr a
x:(r0)-,x0
y:(r4)+,y0
3
rep m0
4
mac x0,y0,a
x:(r0)-,x0
y:(r4)+,y0
5
macr x0,y0,a
(r0)+
6
movep a,x:output
7
jmp fir
This processor has two memory banks called x and y. The code assumes that each in-
put sample can be successively read from a memory location called input, and that
the impulse response is stored in y memory beginning at an address stored in register
r4. Moreover, it assumes that register r0 contains the address of a section of x memory
to use for storing input samples (the delay line), and that this register has been set up
to perform modulo addressing. Modulo addressing means that if it increments or decre-
ments beyond the range of its buffer, then the address wraps around to the other end
of the buffer. Finally, it assumes that the register m0 contains an integer specifying the
number of samples in the impulse response minus one.
The key line (the one that does most of the work) is line 4. It follows a rep
instruction, which causes that one line to be repeatedly executed the number of times
specified by register m0. Line 4 multiplies the data in register x0 by that in y0 and adds
the result to a (the accumulator). Such an operation is called a multiply and accumulate
or mac operation. At the same time, it loads x0 with an input sample from the delay
line, in x memory at a location given by r0. It also loads register y0 with a sample from
the impulse response, stored in y memory at a location given by r4. At the same time,
it increments r0 and decrements r4. This one-line instruction, which carries out several
operations in parallel, is obviously highly tuned to FIR filtering. Indeed, FIR filtering is
a major function of processors of this type.
402
Lee & Varaiya, Signals and Systems
9. FILTERING
x( n)
x( n!1)
x( n!2)
x( n! L+1)
unit
unit
unit
...
delay
delay
delay
Figure 9.14: A delay line.
x( n)
x( n!1)
x( n!2)
x( n! L+1)
unit
unit
unit
...
delay
delay
delay
h(0)
h(1)
h(2)
h( L!1)
y( n)
...
Figure 9.15: A tapped delay line realization of an FIR filter, described as a signal
flow graph.
on end. The right side represents the feedback in the filter, or the recursion that makes it
an IIR filter rather than an FIR filter.
The filter structure in Figure 9.16 can be simplified somewhat by observing that since
the left and right sides are LTI systems, their order can be reversed without changing the
output. Once the order is reversed, the two delay lines can be combined into one, as shown
in Figure 9.17. There are many other structures that can be used to realize IIR filters.
The relative advantages of one structure over another is a fairly deep topic, depending
primarily on an understanding of the effects of finite precision arithmetic. It is beyond the
scope of this text.
9.7
Summary
The output of an LTI system does not contain any frequencies that are not also present
in the input. Because of this property, LTI systems are often called filters. This chapter
considers how filters might be implemented. The convolution sum (for discrete-time sys-
tems) and convolution integral (for continuous-time systems) give descriptions of how to
Lee & Varaiya, Signals and Systems
403
9.7. SUMMARY
x( n)
b(0)
y( n)
unit
unit
delay
b(1)
a(1)
delay
x( n−1)
y( n−1)
unit
unit
delay
b(2)
a(2)
delay
x( n−2)
y( n−2)
...
...
...
...
unit
unit
delay
b( L−1)
a( M−1)
delay
x( n− L+1)
y( n− M+1)
Figure 9.16: A signal flow graph describing an IIR filter. This is called a direct
form 1 filter structure.
calculate an output given an input signal. This description highlights the usefulness of
the response that a system has to an idealized impulse, because the output is given by the
convolution of the input and the impulse response. The frequency response, moreover, is
shown to be Fourier transform of the impulse response, where the Fourier transform is a
generalization of the Fourier series considered in previous chapters. Fourier transforms
are considered in much more depth in the next chapter.
This chapter al