3
DSP Fundamentals
and Implementation
Considerations
The derivation of discrete-time systems is based on the assumption that the signal and
system parameters have infinite precision. However, most digital systems, filters, and
algorithms are implemented on digital hardware with finite wordlength. Therefore DSP
implementation with fixed-point hardware requires special attention because of the
potential quantization and arithmetic errors, as well as the possibility of overflow.
These effects must always be taken into consideration in DSP system design and
implementation.
This chapter presents some fundamental DSP concepts in time domain and practical
considerations for the implementation of digital filters and algorithms on DSP hard-
ware. Sections 3.1 and 3.2 briefly review basic time-domain DSP issues. Section 3.3
introduces probability and random processes, which are useful in analyzing the finite-
precision effects in the latter half of the chapter and adaptive filtering in Chapter 8. The
rigorous treatment of these subjects can be found in other DSP books listed in the
reference. Readers who are familiar with these DSP fundamentals should be able to skip
through some of these sections. However, most notations used throughout the book will
be defined in this chapter.
3.1 Digital Signals and Systems
In this section, we will define some widely used digital signals and simple DSP systems.
The purpose of this section is to provide the necessary background for understanding
the materials presented in subsequent sections and later chapters.
3.1.1 Elementary Digital Signals
There are several ways to describe signals. For example, signals encountered in com-
munications are classified as deterministic or random. Deterministic signals are used
Real-Time Digital Signal Processing. Sen M Kuo, Bob H Lee
Copyright # 2001 John Wiley & Sons Ltd
ISBNs: 0-470-84137-0 (Hardback); 0-470-84534-1 (Electronic)
for testing purposes and for mathematically describing certain phenomena. Random
When the analog sinewave defined in (3.1.3) is connected to the DSP system shown in
Figure 1.1, the digital signal x(n) available for the DSP hardware is the causal sinusoidal
signal
xnA sinVnT f, n 0, 1, ...,I
A sinVnT fun
A sin2pfnT fun,
3:1:5
where T is the sampling period in seconds. This causal sequence can also be expressed
as
xnA sin!n funA sin2pFn fun , 3:1:6
78
DSP FUNDAMENTALS AND IMPLEMENTATION CONSIDERATIONS
where
! VT
V
f
s
3:1:7
is the discrete-time frequency in radians per sample and
F fT
f
f
s
3:1:8
is the normalized frequency to its sampling frequency, f
s
, in cycles per sample.
The fundamental difference between describing the frequency of analog and digital
signals is summarized in Table 3.1. Analog signal sampling implies a mapping of an
infinite range of real-world frequency variable f (or V) into a finite range of discrete-
f cycles per second (Hz) f
F
T
ÀI < f < I
! radians per sample ! 2pF Àp ! p
F cycles per sample F
f
f
s
À
1
2
F
1
2
DIGITAL SIGNALS AND SYSTEMS
79
ynTxn, 3:1:9
where T denotes the computational process for transforming the input signal, x(n), into
the output signal, y(n). A block diagram of the DSP system defined in (3.1.9) is
illustrated in Figure 3.1.
The processing of digital signals can be described in terms of combinations of certain
fundamental operations on signals. These operations include addition (or subtraction),
multiplication, and time shift (or delay). A DSP system consists of the interconnection
of three basic elements ± adders, multipliers, and delay units.
Two signals, x
1
n and x
2
n, can be added as illustrated in Figure 3.2, where
+
Σ
or
Figure 3.2 Block diagram of an adder
x(n) y(n)
a
or
x(n) y(n)
a
Figure 3.3 Block diagram of a multiplier
80
DSP FUNDAMENTALS AND IMPLEMENTATION CONSIDERATIONS
x(n) y(n) = x(n−1)
z
−1
Figure 3.4 Block diagram of a unit delay
is the multiplier's output. The multiply operation of equation (3.1.11) can be imple-
mented by the following C55x code using indirect addressing mode:
amov #alpha, XAR1 ; AR1 points to alpha ()
amov #xn, XAR2 ; AR2 points to x(n)
mpy *AR1, *AR2, AC0 ; AC0 *x(n) y(n)
The sequence {xn} can be shifted (delayed) in time by one sampling period, T,as
illustrated in Figure 3.4. The box labeled z
À1
represents the unit delay, x(n) is the input
signal, and
ynxn À 13:1:12
is the output signal, which is the input signal delayed by one unit (a sampling period).
In fact, the signal xn À 1 is actually the stored signal x(n) one sampling period
(T seconds) before the current time. Therefore the delay unit is very easy to implement
a
(a)
y(n)
Σ
+
+
a
x(n)
x(n−1)
x(n)
x(n−1)
z
−1
z
−1
(b)
Figure 3.5 Block diagrams of DSP systems: (a) direct realization described in (3.1.13), and
(b) simplified implementation given in (3.1.14)
be used to reduce computational requirements. For example, (3.1.13) can be rewritten
as
ynaxnxn À 1: 3:1:14
The block diagram implementation of this difference equation is illustrated in Figure
3.5(b), where only one multiplication is required. This example shows that with careful
design (or optimization), the complexity of the system (or algorithm) can be further
reduced.
The C55x implementation of (3.1.14) can be written as:
amov #alpha, XAR1 ; AR1 points to
amov #temp, XAR2 ; AR2 points to temp
mov *(x1n), AC0 ; AC0 x1(n)
add *(x2n), AC0 ; AC0 x1(n)x2(n)
ynb
0
xnb
1
xn À 1b
2
xn À 2: 3:1:15
The impulse response of the system can be obtained by applying the unit-impulse
sequence dn to the input of the system. The outputs are the impulse response coeffi-
cients computed as follows:
h0y0b
0
Á 1 b
1
Á 0 b
2
Á 0 b
0
h1y1b
0
Á 0 b
1
Á 1 b
2
Á 0 b
1
h2y2b
0
Á 0 b
1
l0
b
l
xn À l
: 3:1:16
Substituting xndn into (3.1.16), the output is the impulse response expressed
as
À 1
hn
LÀ1
l0
b
l
dn À l
b
n
n =0,1,..., L
0 otherwise.
&
3:1:17
Therefore the length of the impulse response is L for the difference equation defined in
(3.1.16). Such a system is called a finite impulse response (FIR) system (or filter). The
impulse response coefficients, b
l
, l 0, 1, ..., L À 1, are called filter coefficients
(weights or taps). The FIR filter coefficients are identical to the impulse response
coefficients. Table 3.2 shows the relationship of the FIR filter impulse response h(n)
and its coefficients b
l
y(n)
Σ
+
+
z
−1
z
−1
b
1
b
0
x(n)
x(n−1)
+
x(n−L+1)
b
L−1
Figure 3.6 Detailed signal-flow diagram of FIR filter
system described by the I/O Equation (3.1.16) is illustrated in Figure 3.6. The string
of z
À1
functions is called a tapped-delay-line, as each z
À1
corresponds to a delay of
one sampling period. The parameter, L, is the order (length) of the FIR filter. The
design and implementation of FIR filters (transversal filters) will be discussed in
Chapter 5.
3.2.1 FIRFilters and Power Estimators
The moving (running) average filter is a simple example of an FIR filter. Averaging is
xnÀxn À L: 3:2:2
Therefore the averaged signal, yn, can be computed recursively as expressed
in (3.2.2). This recursive equation can be realized by using only two additions.
However, we need L 1 memory locations for keeping L 1 signal samples
fxnxn À 1ÁÁÁxn À Lg.
The following C5xx assembly code illustrates the implementation of a moving-
average filter of L 8 based on Equation (3.2.2):
L .set 8 ; Order of filter
xin .usect "indata", 1
xbuffer .usect "indata", L ; Length of buffer
y .usect "outdata", 2,1,1 ; Long-word format
amov #xbufferLÀ1, XAR3 ; AR3 points to endof x buffer
amov #xbufferLÀ2, XAR2 ; AR2 points to next sample
mov dbl(*(y)), AC1 ; AC1 y(nÀ1)in long format
mov *(xin), AC0 ; AC0 x(n)
sub *AR3, AC0 ; AC0 x(n) À x(n±L)
add AC0, #-3, AC1 ; AC1 y(n±1) 1/L[x(n)±x(n±L)]
mov AC1, dbl(*(y)) ; y(n) AC1
rpt # (LÀ1) ; Update the tapped-delay-line
mov *AR2À, *AR3À ; X(n±1) x(n)
mov *(xin), AC0 ; Update the newest sample x(n)
mov AC0, *AR3 ; X(n) input xin
The strength of a digital signal may be expressed in terms of peak value, energy, and
power. The peak value of deterministic signals is the maximum absolute value of the
signal. That is,
M
x
max
n
fjxnjg: 3:2:3
L
LÀ1
n0
jxnj
2
: 3:2:5
If xn is a periodic signal, we have
xnxn kL, 3:2:6
where k is an integer and L is the period in samples. Any one period of L samples
completely defines a periodic signal. From Figure 3.7, the power of xn can be
computed by
P
x
1
L
n
lnÀL1
jxlj
2
1
L
LÀ1
l0
jxn À lj
2
x
n À 1
1
L
x
2
nÀx
2
n À L: 3:2:9
86
DSP FUNDAMENTALS AND IMPLEMENTATION CONSIDERATIONS
To further simplify the algorithm, we assume L is large enough so that
x
2
n À L%
^
P
x
n À 1 from a statistical point of view. Thus Equation (3.2.9) can be
further simplified to
^
P
x
n% 1À
1
L
^
P
x
for non-stationary signals for better results. In many real-time applications, the square
of signal x
2
n used in (3.2.10) and (3.2.11a) can be replaced with its absolute value jxnj
in order to reduce further computation. This efficient power estimator will be further
analyzed in Chapter 4 using the z-transform.
3.2.2 Response of Linear Systems
As discussed in Section 3.1.3, a digital system can be completely characterized by its
impulse response hn. Consider a digital system illustrated in Figure 3.8, where xn is
the input signal and yn is the output signal. If the impulse response of the system is
hn, the output of the system can be expressed as
ynxnÃhn
I
kÀI
xkhn À k
I
kÀI
hkxn À k, 3:2:12
x(n) y(n) = x(n)∗h(n)
h(n)
Figure 3.8 A simple linear system expressed in time domain
INTRODUCTION TO DIGITAL FILTERS
87
where * denotes the linear convolution operation and the operation defined in (3.2.12) is
called the convolution sum. The input signal, xn, is convoluted with the impulse
response, hn, in order to yield the output, yn. We will discuss the computation of
linear convolution in detail in Chapter 5.
As shown in (3.2.12), the I/O description of a DSP system consists of mathematical
A digital filter can be classified as either an FIR filter or an infinite impulse response
(IIR) filter, depending on whether or not the impulse response of the filter is of
finite or infinite duration. Consider the I/O difference equation of the digital system
expressed as
ynbxnÀayn À 1, 3:2:16
where each output signal yn is dependent on the current input signal xn and the
previous output signal yn À 1. Assuming that the system is causal, i.e., yn0 for
n < 0 and let xndn. The output signals yn are computed as
88
DSP FUNDAMENTALS AND IMPLEMENTATION CONSIDERATIONS
y0bx0ÀayÀ1b,
y1bx1Àay0Àay0Àab,
y2bx2Àay1Àay1a
2
b,
...
In general, we have
hnynÀ1
n
a
n
b, n 0, 1, 2, ...,I: 3:2:17
This system has infinite impulse response hn if the coefficients a and b are non-zero. This
system is called an IIR system (or filter). In theory, we can calculate an IIR filter output
yn using either the convolution equation (3.2.14) or the I/O difference equation (3.2.16).
However, it is not computationally feasible using (3.2.14) for the impulse response hn
given in (3.2.17), because we cannot deal with an infinite number of impulse response
coefficients. Therefore we must use an I/O difference equation such as the one defined in
(3.2.16) for computing the IIR filter output in practical applications.
The I/O equation of the IIR system given in (3.2.16) can be generalized with the
, m 1, 2, ..., Mg. Since the
outputs are fed back and combined with the weighted inputs, this system is an example
of the general class of feedback systems. Note that when all a
m
are zero, Equation
(3.2.18) is identical to (3.1.16). Therefore an FIR filter is a special case of an IIR filter
without feedback coefficients. An FIR filter is also called a non-recursive filter.
The difference equation of IIR filters given in (3.2.18) can be implemented using the
MATLAB function filter as follows:
yn filter(b, a, xn);
where the vector b contains feedforward coefficients fb
l
, l 0, 1, ..., L À 1g and the
vector a contains feedback coefficient fa
m
, m 1, 2, ..., Mg. The signal vectors, xn
and yn, are the input and output buffers of the system. The FIR filter defined in (3.1.16)
can be implemented using MATLAB as
yn filter(b, 1, xn);
Assuming that L is large enough so that the oldest sample xn À L can be approxi-
mated using its average, yn À 1. The moving-average filter defined in (3.2.2) can be
simplified as
yn 1À
1
L
yn À 1
1
L
xn1 À ayn À 1axn, 3:2:19
fx Xg. Some properties of F(X) are summarized as follows:
FÀI 0, 3:3:4a
FI 1, 3:3:4b
0 FX 1, 3:3:4c
90
DSP FUNDAMENTALS AND IMPLEMENTATION CONSIDERATIONS
FX
1
FX
2
if X
1
X
2
, 3:3:4d
PX
1
< x X
2
FX
2
ÀFX
1
: 3:3:4e
The probability density function of a random variable x is defined as
fX
dFX
dX
3:3:5
if the derivative exists. Some properties of f(X) are summarized as follows:
0, x < X
1
or x > X
2
a, X
1
x X
2
,
&
which is uniformly distributed between X
1
and X
2
. The constant value a can be
computed by using (3.3.6b). That is,
I
ÀI
fXdX
X
2
X
1
a Á dX aX
2
À X
1
1:
−X
1
Figure 3.9 A uniform density function
fX
1
X
2
À X
1
, X
1
X X
2
0, otherwise.
@
3:3:7
This uniform density function will be used to analyze quantization noise in Section 3.4.
If x is a discrete random variable that can take on any one of the discrete values
X
i
, i 1, 2, ... as the result of an experiment, we define
p
i
Px X
i
: 3:3:8
3.3.2 Operations on Random Variables
We can use certain statistics associated with random variables. These statistics
areoften more meaningful from a physical viewpoint than the probability density
function. For example, the mean and the variance are used to reveal sufficient
follows:
92
DSP FUNDAMENTALS AND IMPLEMENTATION CONSIDERATIONS
Table 3.3 Probability of rolling a die
X
i
123456
p
i
1/6 1/6 1/6 1/6 1/6 1/6
The mean of outcomes can be computed as
m
x
6
i1
p
i
X
i
1
6
1 2 3 4 5 63:5:
The variance of x, which is a measure of the spread about the mean, is defined as
s
2
x
Ex À m
. The mean of the squared
deviation indicates the average dispersion of the distribution about the mean m
x
. The
positive square root s
x
of variance is called the standard deviation of x. The MATLAB
function std calculates standard deviation. The statement
s std(x);
computes the standard deviation of the elements in the vector x.
The variance defined in (3.3.10) can be expressed as
s
2
x
Ex À m
x
2
Ex
2
À 2xm
x
m
2
x
Ex
2
À2m
x
Exm
y
a
2
s
2
x
. It can also be shown (see exercise problem) that P
x
m
2
x
s
2
x
if m
x
T 0.
Consider a uniform density function as given in (3.3.7). The mean of the function can
be computed by
INTRODUCTION TO RANDOM VARIABLES
93
m
x
Ex
I
ÀI
XfXdX
1
X
X
2
fXdX À m
2
x
1
X
2
À X
1
X
2
X
1
X
2
dX À m
2
x
1
X
2
À X
1
:
X
3
2
À X
1
2
12
:
3:3:14
In general, if x is a uniformly distributed random variable in the interval (ÀD, D), the
mean is 0 (m
x
0) and the variance is s
2
x
D
2
=3.
Example 3.3: The MATLAB function rand generates pseudo-random numbers
uniformly distributed in the interval (0, 1). From Equation (3.3.13), the mean of
the generated pseudo-random numbers is 0.5. From (3.3.14), the variance is
computed as 1/12. To generate zero-mean random numbers, we subtract 0.5
from every generated number. The numbers are now distributed in the interval
(À0.5, 0.5). To make these pseudo-random numbers with unit-variance, i.e.,
s
2
x
D
2
=3 1, the generated numbers must be equally distributed in the interval
(À
random variable y is defined as
y x
1
x
2
ÁÁÁx
N
N
i1
x
i
: 3:3:17
94
DSP FUNDAMENTALS AND IMPLEMENTATION CONSIDERATIONS
m
y
y
f(y)
1/s
y
√2p
Figure 3.10 Probability density function of Gaussian random variable
The probability density function f(Y ) becomes a Gaussian (normal) distribution func-
tion (normal curve) as N 3I. That is,
fY
1
s
y
where m
y
N
i1
m
i
and s
y
N
i1
s
2
i
q
. A graphical representation of the probability
density function defined in (3.3.18) is illustrated in Figure 3.10.
The central limit theorem defined in (3.3.17) is useful in generating a Gaussian
random variable from a uniformly distributed random variable using N ! 12. The
Gaussian random variable is frequently used in communication theory. The MATLAB
function randn generates pseudo-random numbers normally distributed with mean 0
and variance 1.
3.4 Fixed-Point Representation and Arithmetic
The basic element in digital hardware is the two-state (binary) device that contains one
bit of information. A register (or memory unit) containing B bits of information is called
a B-bit word. There are several different methods of representing numbers and carrying