INTRODUCTION TO

DIGITAL SIGNAL PROCESSING

"This page is Intentionally Left Blank"

INTRODUCTION
TO DIGITAL
SIGNAL PROCESSING

RomanKuc
Department of Electrical Engineering
Yate University

SSP BS Publications
4-4-309, Giriraj Lane, Sultan Bazar,
Hyderabad - 500 095 - AP.
Phone: 040-23445677,23445688
e-mail: ci.: . [email protected]
www.bspublications.net

Reprint 2008
All rights reserved.

No part of this book or parts thereof may be

reproduced, stored in a retrieval system or
transmitted in any language or by any means, ,
electronic, mechanical, photocopying, recording ,
or otherwise without the prior written
permission of the publishers.
I

Copyright 1982, by Author

"This Reprint of Introduction to Digital Signal
Processing, originally published by McGraw Hill, Inc.
1982, is published by agreement with the author".

Published by :

SSP BS
=
__

Publications

4-4-309, Giriraj Lane, Sultan Bazar,

Hyderabad - 500 095 - A. P.
Phone: 040-23445677,23445688
Fax No: 040-23445611

e-mail: [email protected]
www.bspublications.net

Printed at
Adithya Art Printers

Ilyderabad

ISBN

81-7800-168-3

PREFACE

This book is intended to be used in the first course covering digital signal
processing and filter design, typically offered at the senior or first-year
graduate level in electrical engineering. The course is also appropriate for a
graduate course in departments other than electrical engineering, such as
geophysics and mechanical engineering, in which the analysis of discrete"time
data is performed. It is assumed that the student has had a course covering
Fourier series and LaPlace transforms on the level of the first linear circuits or
control systems course. This text also includes projects that require students to
write computer programs to accomplish signal processing projects. The student
should bO'familiar with some programming language 'such as FORTRAN,
BASIC, PASCAL, or C.
This book approaches digital signal processing and filter design in a novel
way, by presenting the relevant theory and then having the student apply it by
implementing signal processing routines on a computer. This mixture of theory
and application has worked successfully for the past six years' in teaching this
course at Yale University. With this approach, the students receive a deeper
and intuitive understanding of the theory, its applications and its limitations.
The text can accommodate a wide variety of courses. Currently, the
course on digital filters at the undergraduate level are taught primarily as a
theory course, with homework problems and exams to determine the course
grade. This book can be used directly in this type of course. However, the
course becomes much more interesting and enjoyable when students apply the
theory to write programs to perform signal processing tasks. Each chapter
describes programming techniques that implement the theory and includes
projects that illustrate different applicatio'!l:>. Digital filter routines are written
from the very beginning of the course, immediately after the time-domain
description of discrete-time systems is presented in Chapter 2. The programs
can be written on any computer from personal to mainframe, requiring only a
standard terminal for displaying the results in graphical form and a printer for
v

vi

PREFACE

generating the hard copy. The programs to be written by the student are not
difficult. Each routine performs only one primitive task in the signal processing
procedure. For example, the discrete Fourier transform subprogram can be
accomplished in only nine lines of code. But once the student has given
sufficient thought to the problem to write these nine lines, the magic associated
with this procedure disappears, being replaced by understanding. To minimize
the programming effort of the student, the routines are used as building blocks
in later projects. Constructing more sophisticated programs from these simpler
building blocks aids the student in conceptualizing projects later in the course.
By the end of the course, the student has a package of signal processing
programs that can be used for a variety of applications.
The graphic outputs for the projects are displayed on the terminal and on
a printer, being generated by programs provided in the Appendix. These plots,
although having crude resolution, are adequate for illustrating the ideas and
also have the advantage of not requiring special plotting hardware. Students
having higher-resolution terminal graphic capabilities have readily written their
own graphic routines using the programs in the Appendix as guides.
This text differs from others in the area of digital filters and signal
processing in the important respect that processing signals and designing filters
on a digital computer, using simple programs composed by the students, are
an integral part of the course. For this reason, more realistic problems can be
assigned and discussed than would be in a course without computers, for whichthese problems would be mathematically untenable or at least tedious. The
organization of the material presented in the text allows the student to start
writing meaningful filtering programs from the beginning of the course.
Projects included in the text illustrate both digital signal processing and digital
filter design concepts.
The book is organized in the following manner. Chapter 1 defines the
basic components of a digital filter and provides motivation by describing
several common applications of digital filters and signal processing. The
material in Chapters 2 through 5 covers the fundamentals of time, frequency,
and z-plane analysis. If some of this material has been introduced in previous
courses, the respective chapters can be treated as a review. Several common
structures that are employed to implement digital filters are described in
Chapter 6.
The synthesis of digital filters to meet a desired specification is presented
in Chapters 7 through 9. To compare the different filter design procedures, the
projects are structured to allow each student to be assigned a personal
specification to be satisfied. In Chapter 7, the digital filter design is accomplished by interactively and intuitively placing poles and zeros in the z plane.
This procedure serves three purposes. First, it develops an understanding of
the interaction of the singularities in generating a filter magnitude response or
system spectrum. Second, it provides an appreciation for the conventional
filter design procedures presented in Chapters 8 and 9. Finally, this intuitive
procedure indicates how to correct the locations of the poles and zeros

PREFACE

vii

generated by conventional design procedures. The computed positions of these

singularities are slightlY, but significantly, displaced from the desired positions
by the round-off errors that are always present in any computation. These
round-off effects become obvious in the filter design projects. The conventional
procedures for designing infinite-impulse response (IIR) digital filters are
described in Chapter 8, and those for finite-impulse response (FIR) digital
filters in Chapter 9.
Two advanced topics are covered in Chapters 10 and 11. Chapter 10
describes the effects of finite precision encountered in practical hardware
implementation of digital filters. Some elementary results from probability
theory are presented to describe noise effects. Chapter 11 provides an
overview of inverse filters. This topic was chosen because it includes and
integrates many of the concepts introduced in the previous chapters to
illustrate an important practical application.
In addition to the Problems, a set of Computer Projects is included, the
programs for which should be written by the student. Instructions for program
development are included as an integral part of the text and Project
description. The ability to manipulate sequences by a computer program is
essential for appreciating the concepts involved for signal processing, which is
important when the studen~ IS faced with novel problems in his or her own
application.
Projects in the analysis part of the book illustrate the practical applications of digital interpolation procedures and the simulation of a frequency-shift
keying (FSK) digital communication channel. In the filter synthesis part of the
book, the Projects include designing digital filters to satisfy a practical
specification. In the chapters describing finite-precision effects and the inverse
filter, the Projects illustrate programming techniques for simulating physical
hardware and systems.
One main benefit of the Projects is to relate the analytic results to those
produced by the computer. At times, the analysis is simple and should then be
used to predict and verify the computer results. At other times, the analytic
procedure is very complicated mathematically, in which case the computer
results should be obtained first to indicate what type of results should be
expected. In this latter case, the computer results indicate the direction in
which the analysis should proceed. For example, the computer results may
indicate that a complex-valued operation in a particular analytic step results in
a real-valued sequence. This knowledge may serve to drive the analytic steps.
In practice, it is often found that computer simulations demonstrate an effect
that can be proven analytically.
An insider's knowledge of the signal processing algorithms is attained by
having the students write their own programs. This knowledge then allow~ the
student comfortably to make the intellectual move to more sophisticated signal
processing and filter design programs and packages, such as those distributed
by the IEEE or SIG or offered by many vendors.
This text has the flexibility of being used in a wide variety of digital signal

viii

PREFACE

processing courses at the undergraduate and graduate levels. For a conventional undergraduate theory course in digital filters, the first nine chapters of
this book supply the necessary topics in discrete-time signals and systems. In
some universities, such a theory course is followed by a digital filter design
laboratory. The Projects in this book can be used as the exercises in such a
course. At the graduate level, a faster pace allows topics from all chapters to
be covered. Courses stressing the signal processing aspects should emphasize
the analytic methods and the discrete Fourier transform and de-emphasize IIR
filter design and quantization effects, while courses stressing filter implementation should stress the latter and reduce the emphasis on the former.
To allow adequate space for this treatment of this material, the analytic
techniques have been restricted to those that are absolutely necessary for
understanding the filter synthesis and signal processing procedures. Topics that
are normally covered in a senior/first-year graduate course on digital signal
processing, but that are not considered in this text, include the solution of the
difference equations with the one-sided z-transform, the numerical methods
for optimizing digital filter designs, and the area of multidimensional signal
processing. These topics are better treated in a second course in which the
students have the familiarity of the material contained in this book.
The author is indebted to the students at Yale University who have taken
this course and provided many useful comments regarding the topics and the
exercises in this book. Among these are Ata Arjomand, Billur Barshan, Leo
Costello. An;lOld Englander, Kannan Parthasarathy, Rob Roy, and David
Starer. The author would like to thank the following reviewers for their
comments and-suggestions, J. I. Aunon, Purdue University; Edward Delp,
Purdue University; Simon Haykin; McMaster University; M. Kaveh, University of Minnesota; and Paul Prucnal, Columbia University.
The author hopes that the instructors teaching digital signal processing as
a combination of theory and computer exercises, as described in this text, have
an enjoyable learning experience. This learning comes about from the
questions posed by the students, who inevitably find some novel approach or
interpretation to filter design on the computer. This has certainly been the case
every time the author has taught this course.
Roman Kuc

CONTENTS

Preface
1 Introduction
1-1
1-2
1-3
1-4
1-5
1-6
1-7
1-8

What is a Digital Filter?

Anatomy of a Digital Filter
Frequency Domain Description of Signals and Systems
Some Typical Applications of Digital Filters
Replacing Analog Filters with Digital Filters
Overview of the Course
Description of the Computer Projects
Summary

v
1
1
3
5
7
12
13
17
18

2 Discrete-time Description of Signals and

Systems
2-1
2-2
2-3
2-4
2-5
2-6
2-7
2-8
2-9
2-10
2-11

Introduction
Discrete-Time Sequences
Superposition Principle for Linear Systems
Unit-sample Response Sequence
Time-Invariant Systems
Stability Criterion for Discrete-time Systems
Causality Criterion for Discrete-time systems
Linear Constant-Coefficient Difference Equations
Suggested Programming Style
Writing a Digital Filter Program
Summary

3 Fourier Transform of Discrete-time Signals

3-1
3-2
3-3

Introduction
Definition of the Fourier Transform
Important Properties of the Fourier Transform

20
20
20
26
27
30
38
40

41
44

45
55
65
65
65
69
ix

CONTENTS

3-4
3-5
3-6
3-7
3-8
3-9
3-10
3-11
3-12

4
4-1
4-2
4-3
4-4
4-5
4-6
4-7
4-8
4-9
4-10
4-11
4-12
4-13
4-14

Properties of the Fourier Transform for Real-valued Sequences

Program to Evaluate the Fourier Transform by Computer
Use of the Fourier Transform in Signal Processing
Fourier Transform of Special Sequences
The Inverse Fourier Transform
Fourier Transform of the Product of Two Discrete-time
Sequences
Sampling a Continuous Function to Generate a Sequence
Reconstruction of Continuous-time Signals from Discrete-time
Sequences
Summary

The Discrete Fourier Transform

Introduction
The Definition of the Discrete Fourier Transform (DFT)
Computing the Discrete Fourier Transform from the Discretetime Sequence
Properties of the DFT
Circular Convolution
Performing a Linear Convolution with the DFT
Computations for Evaluating the DFT
Programming the Discrete Fourier Transform
Increasing the Computational Speed of the DFT
Intuitive Explanation for the Decimation-in-time FFT
Algorithm
Analytic Derivation of the Decimation-in-time FFT Algorithm
Some General Observations about the FFT
Other Fast Realizations of the DFT
Summary

5 The z-transform
5-1
5-2
5-3
5-4
5-5
5-6
5-7
5-8
5-9
5-10
5-11
5-12
5-13
5-14

Introduction
The Definition of the z-transform
Properties of the z-transform
The System Function of a Digital Filter
Combining Filter Sections to Form More Complex Filters
Digital Filter Implementation from the System Function
The Complex z-plane
The Region of Convergence in the z-plane
Determining the Filter Coefficients from the Singularity
Locations
Geometric Evaluation of the z-transform in the z-plane
. Relationship Between the Fourier Transform and the
z-transform
The z-transform of Symmetric Sequences
The Inverse z-transform
Summary

72
80

82
84
92
94
96
105
110
118
118
119
121
128
130
132
133
135
136
142
146
148
150
150
157
157
158
160
162
164
167
169
i74
178
183
184
188
196
203

CONTENTS

6 Digital Filter Structures

6-1
6-2
6-3
6-4
6-5
6-6
6-7
6-8
6-9

Introduction
System Describing Equations
Filter Catagories
The Direct Form I and II Structures
Cascade Combination of Second-order Sections
Parallel Combination of Second-order Sections
Linear-phase FIR Filter Structures
Frequency-sampling Structure for the FIR Filter
Summary

7 From Analysis to Synthesis

7-1
7-2
7-3
7-4
7-5
7-6
7-7
7-8
7-9
7-10
7-11

Introduction
Interrelationships of Analytic Methods
Practical Magnitude Response Specifications
Log-magnitude Response Curves
Programs to Compute the Log-magnitude Response of Poles
and Zeros
Phase Response Considerations
Steps in Performing a Filter Design
Choice of Filter Type
Interactive Filter Design by Intuitive Pole/Zero
Placement
Writing an Interactive Design Program
Summary

8 Infinite Impulse Response Filter Design

Techniques
8-1
8-2
8-3
8-4
8-5
8-6
8-7

Introduction
Analog Filter System Function and Frequency Response
Analog Lowpass Filter Design Techniques
Methods to Convert Analog Filters into Digital Filters
Frequency Transformations for Converting Lowpass Filters
into other Types
All-pass Filters for Phase Response Compensation
Summary

9 Finite Impulse Response Filter Design

Techniques
9-1
9-2
9-3
9-4
9-5

Introduction
Designing FIR Filters with the Windowing Method
The DFT Method for Approximating the Desired Unit-Sample
Response
Designing FIR Filters by Combining the DFT and Window
Methods
Designing FIR Filters with the Frequency-sampling Method

xi

208
208
209
209
210
214
218
220
221
227
230
230
231
236
238
239
243
246
246
247
256
258

262
262
263
263
280
295
310
312

318
318
319
341
348

361

xii

CONTENTS

9-6
9-7

Use of the FFT Algorithm to Perform FIR Filtering

Summary

10 Finite-precision Effects
10-1
10-2
10-3
10-4
10-5
10-6
10-7
10-8

Introduction
Finite-precision Number Representation within the Computer
Quantization Error in Analog-to-digital Conversion
Expressing Coefficients in Finite Precision
Performing Arithmetic in Finite-Precision Number Systems
Effects of Finite-Precision Arithmetic on Digital Filters
Programs to Simulate Quantization Effects
Summary

11 Inverse Filtering
11-1
11-2
11-3
11-4
11-5
11-6
11-7
11-8
11-9

Introduction
Applications of Inverse Filters
Minimum Phase Systems
Formulating the Problem for Applying an Inverse Filter
FIR Filter Approximations to the Inverse Filter
Discrete Hilbert Transform Relationship for Minimum-Phase
Systems
Designing Inverse Filters with the Discrete Hilbert Transform
The Effects of Noise on Inverse Filters
Summary

379
384
390
390
391
395
399
404
409

411
416
422
422
423
426
429
432
433
436
441
448

Appendix A

455

Appendix B

469

Index

471

CHAPTER

1
INTRODUCTION

In this chapter, we present the concept of a digital filter and define the adder,
multiplier and delay elements that are its components. Some typical applications of digital filters for signal processing, system simulation and spectral
analysis are described. An overview of the topics covered in the book
concludes the chapter.

1.1 WHAT IS A DIGITAL FILTER?

Most broadly defined, a digital filter is a numerical procedure, or algorithm,
that transforms a given sequence of numbers into a second sequence that has
some more desirable properties, such as less noise or distortion. As shown in
Fig. 1.1, the entire given, or input, sequence will be denoted in this book by
{x(n}}, and the entire output sequence by (y(n)}, where n is an index.
Typically, the set of index values consists of consecutive integers, which in
some cases may take values from minus infinity to plus infinity.
The desired features in the output sequence depend on the application.
For example, if the input sequence is generated by a sensing device, such as a
microphone, the digital filter may attempt to produce an output sequence
having less background noise or/interference. In radar applications, digital
filters are used to improve the5detection of airplanes. In speech processing,

INTRODUCTION TO DIGITAL SIGNAL PROCESSING

1""____
. .I----o...,uIPul{y(n))sequence
DIGITAL

Input seQueo-nce
__
_
{*)} 0

FILTER

FIGURE I-I
A digital filter transforms the input
sequence {x(n)} into the output
sequence (y(n)}.

digital filters have been employed to reduce the redundancy in the speech
signal so as to allow more efficient transmission, and for speech recognition.
Input sequences can be generated in several ways, as shown in Fig. 1.2.
One common method is to sample a continuous-time signal at a set of
equally-spaced time intervals. If the continuous-time signal is denoted by
x(t), then the values of the discrete-time sequence are denoted as
x(nT.) = x(t)lt=nT,

(1.1)

where T. is the sampling period. In image processing applications, scanning a

line of an image produces an intensity function of position across the image
plane, denoted by i(p). This function can be sampled to produce a discreteposition sequence with values given by
i(nD.) = i(P )lp=nD.

(1.2)

where D. is the sampling distance interval.

One common application of digital filters is to simulate a physical system.
The input sequence is then usually a sequence of numbers generated by the
computer and represents the excitation signal {e(n)}. In this case, {e(n)}
represents instantaneous values of some physical parameter, such as VOltage,
torque, temperature, etc. This sequence is then processed to determine the

CONTINUOUS-TIME SIGNAL

DISCRETE-TIME SEQUENCE
{x(nT.)}

fIII,

IMAGE

III I JlI ,.

DISCRETE-TIME SEQUENCE
{i(nD.)}

-+#-iIR-IllII--IJII--t-Scan line

....IUI.....rr.......

FIGURE 1-%
Deriving one-dimensional discrete-time sequences from analog signals.

nD.

nT.

INTRODUCTION

response of the system, represented by a sequence {r(n)}. This application is

discussed in more detail later in this chapter.
To allow these different sequence types to be treated with the same
notation, we will generalize the element index to be an integer, denoted by n.
That is, we will write x(n), rather than x(nT.), and i(n), rather that i(nDs).
In addition to sampling, the input sequence could have come in a
discrete-time form initially, such as daily sales receipts, {m(n)}, or weekly
measurements of personal weight, {w(n)}. In these cases, we may be
interested in short-term fluctuations or long-term trends.
The above data sequences are indexed in terms of temporal or spatial
increments. To simplify the discussion of the various types of input data
sequences that can exist, we will refer to all of them as discrete-time sequences.
This is done with the understanding that the particular sequence under
discussion may, in practice, be related to a physical parameter other than time.
Although the theory presented in this text is general, the implementation
of a digital filter depends on the application and the user's environment. In
education and research, a digital filter is typically implemented as a program
on a general-purpose computer. The types of computers that can be used vary
widely, from personal computers, through the larger minicomputers, to the
large time-shared mainframes. In commercial instrumentation and industrial
applications, the digital filter program is commonly implemented with a
mierocomputer that may also be used for control and monitoring purposes as
well. For high-speed or large-volume applications, such as use in automobiles
for controlling the engine operation, the digital filter may consist of specialpurpose integrated circuit chips. These circuits perform the computation and
storage functions required for the digital filter operation. These functions,
which are common to all digital filter applications, are described in the next
section.

1.2 ANATOMY OF A DIGITAL FILTER

A digital filter consists of the interconnection of three simple elements: adders,
multipliers and delays. These are shown in Fig. 1.3. The adder and multiplier
are conceptually simple components that are readily implemented in the
arithmetic logic unit of the computer. Delays are components that allow access
to future and past values in the sequence. Open-headed arrows indicate
direction of information flow, while the larger c1o~ed-headed arrows indicate
multipliers. This convention will be useful for drawing complicated digital
filters in the later chapters.
Delays come in two basic "flavors": positive and negative. A positive
delay, or simply delay, is implemented by a memory register that stores the
current value of a sequence for one sample interval, thus making it available
for future calculations. A positive delay is conventionally indicated by a box
denoted by Z-l. The significance of this notation will become clear when we
analyze digital filters with the z-transform. A negative delay, or advance, is

INTRODUCTION TO PlGrrAL SIGNAL PROCESSING

ADDER

---;-cp"
+ "'--_0

xl(n) ...

y(n)

=xl(n) + x2(n)

x2(n)-

MULTIPLIER

POSITIVE AND
NEGATIVE DELAYS

x(n)

K"'---__

----iV-

_0

y(n) = ax(n)

POSITIVE DELAY - "DELAY"

--.".B. .---o
8

x(n) ....

y(n) = x(n-I)

NEGATIVE DELAY - "ADVANCE"

x(n) 0

y(n) = x(n+ I)

nGURE 1-3
Elements that are interconnected to implement a digital filter.

used to look ahead to the next value in the sequence, and is indicated by a box

denoted by z. Advances are typically used for applications, such as image

processing, in which the entire data sequence to be processed is available at
the start of processing, so that the advance serves to access the next data
sample in the sequence. The availability of the advance will simplify the
analysis of digital filters. However, the advance cannot always be implemented
in some appliCations. For example, when the data sequence is obtained by
sampling a function of time, each new sample is usually processed as it is
acquired. In this case, advances are not allowed, since we cannot gain access to
future data values.
A digital filter, design involves selecting and interconnecting a finite
number of these elements and determining the multiplier coefficient values.
Examples of interconnecting these elements to form two simple digital filters
are considered in the following examples.
EumpIe 1.1. Digital filter used as a tbree-sample averager. Consider the
relationship between the values of the output sequence at time n, denoted by
y(n), and the values of the input sequence, x(n -1), x(n) and x(n + 1), given by

y(n) = ix(n

+ 1) + lx(n) + ix(n -1)

As shown in Fig. 1.4, the current output y(n) is equal to the average of the next,
current, and previous input values. The advance serves to access the next value of
the sequence, while the delay stores the previous value.

Any past or future input data values can be accessed by including the
appropriate number of delay or advance elements. Data values can be
combined, however, only through an ad~er.

INTRODUcnON

x(n)--+------..;-c

y(n) = ih [x(n+ I) + x(n)+ x(n-I)]

FIGURE 1....
Nonrecursive digital filter that acts as a three-sample averager.

The term order will be used to denote the minimum number of delays
and advances that are necessary to implement a digital filter. For example, the
3-sample averager is a second-order filter. When the filter output is a function
of only the input sequence, as in the 3-sample averager, it is called a
nonrecursive filter. When the output is also a function of the previous output
values, it is called a recursive filter. Example 1.2 illustrates a recursive filter.
Example 1.2. First-order recursive filter. Consider the relationship between the
input and output sequence values given by
y(n) = ay(n -1)

+ x(n)

where a is a constant. As shown in Fig. 1.5, the current value of the output is
equal to the sum of the input and a times the value of the previo,lls output. Since
only one delay is necessary, this is a first-order filter.

Example 1.2 shows how feedback is implemented in a digital filter. Whereas

both positive and negative delays can be applied to the input sequence, only
positive delays can be applied to the output sequence, since we cannot gain
access to output values that have not yet been computed.

1.3 FREQUENCY DOMAIN

DESCRIPI10N OF SIGNALS AND
SYSTEMS
To discuss some typical applications of digital filters in conventional terms, we
need to define the frequency domain description of signals and systems. Let us

"')If-::-(--a--$~Z_I--
'----l,\r--<<J
FIGURE 15
A first-order recursive digital filter.

y(n) =ay(n-I) + x(n)

(;

INTRODUcnON TO DIGITAL SIGNAL PROCESSING

(x(n)} =

{COS

(wn)}

".

"

nGURE 1-6
The cosine sequence {cos(wn)} for w = :tr/l2.

. DISCRETE-TIME SEQUENCES

MAGNITIJDE SPECrRA

~--~-----:!:~w

(X2(n))

... 1'

I II I I I I

JI

1...

)n
oL _ _ _

---L~

_ _~.;. W

I
nGURE 17
Discrete-time sequences and their magnitude spectra. The sequence {x)(n)} = {x\(n)} + {x2(n)}.

7
consider the discrete-time sinusoidal sequence {x(n)} shown in Fig. 1.6, and
given by
INTRODUcrlON

{x(n)}

= {cos(wn)}

(1.3)

where w is the frequency of the sequence, having units of radians. More

interesting signals can be considered to be composed of a sum of such
sinusoidal sequences having different frequencies. One example is the Fourier
series expansion of a periodic waveform. As described in detail in Chapter 3,
the magnitude spectrum of a discrete-time sequence, denoted by IX(eIW)1
indicates the magnitude of the various sinusoidal components. Examples of
three sequences and their magnitude spectra are shown in Fig. 1.7. The
sequence {xl(n)} contains mostly low-frequency components, {x2(n)} contains
mostly high-frequency components, and {x3(n)} = {xl(n)} + {x2(n)} contains
both.
One common application of digital filters is to separate signals on the
basis of their spectral content. For example, in Fig. 1.8, we start with the
sequence {x3(n)}. To retrieve the low-frequency signal {xl(n)}, we pass
{x3(n)} through a filter whose magnitude response, denoted by IHLP(eiW)I,
passes the low-frequency components of the input, while eliminating the
high-frequency ones. Such a filter is called a lowpass filter. To retrieve the
high-frequency signal {x2(n)}, we pass {x3(n)} through a highpass filter, whose
response is shown as IHHP(eiw)l. Two other important filters are the bandpass
filter, which passes the input components that lie only within a desired band of
frequencies, and the bandreject filter, which eliminates the input components
that lie within a specified band of frequencies.

1.4 SOME TYPICAL APPLICATIONS OF

DIGITAL FILTERS
Studying the field of digital filtering provides a fundamental understanding of
processing of discrete-time data. This is essential for manipulating data by
digital computer and for applying the signal processing procedures to such
applications as robotics, intelligent sensors and seismic signal processing.
Digital filter teChniques are also important for modeling linear systems, such as
speech production models and transmission channels. Often, such modeling
provides insights into the underlying principles of the problem, from which
creative solutions can stem. We now consider some typical applications.
SPEECH PROCESSING. Digital filters are commonly used for speech recognition systems. Let us consider the sequences {xl(n)} and {x2(n)}, shown in Fig.

1.9, which represent samples taken from two speech signals. The first is the lal
sound, as in "father", and the second is the IiI sound, as in "see". These two
periodic signals are most easily differentiated by computing the magnitude
spectrum of one period. We note that the two spectra, IX1(eiw)1 and IX2(e IW)I,
consist of peaks that correspond to the resonances of the vocal tract shape that

INTRODUcnON TO DIGITAL SIGNAL PROCESSING

{X3(n) }

.,1"1,, III

...

1...

IHuo(ei"') I

{X3(n)}~---J

LOWPASS

L{yu.(n)}

I . . ._FI_I:_:rE_R_-,1 =

{x.(n)}

'I

{X3(n)}~ ~IGHPASS L{y~n)}

FILTER

r--

{x~n)}

I..

F:GURE 1-8
Recovering sequences by filtering. The sequence {Xt(n)} is recovered from {x3(n)} by lowpass
filtering. and {x2(n)} is recovered by bighpass filtering.

was used to generate the sound. In speech recognition systems, various speeCh
sounds can be differentiated on the basis of the frequency locations of these
peaks.
COMPUTER SIMULATION OF A PHYSICAL SYSTEM. When the physical
system, such as a robot, is very complicated or expensive to build, often it is
first simulated on a digital computer in its digital filter form. This simulation is
employed to determine the performance and the sensitivities of the analog

INTRODUcnON

(X.(II)}

-~~O------------------------------~nro

Pvr\A,,,.

(X,(II)}

~~O------------------------------~nro

1'

FIGURE
Speech signals and their spectra. The sequence {x1(n)} consists of 512 samples of the speech sound
lal and {x2(n)} of IiI. The sample values are connected by straight lines. The corresponding
spectra IX1(ejO1 and IX2(ei"')I, given in logarithmic units of dB, exhibit peaks that can be used to
differentiate these sounds.

10

INTRODUCTION TO DIGITAL SIGNAL PROCESSING

system before it is actually built. The frequency-dependent transformations,

such as those produced by the inertia and compliance of a robot arm, can be
mimicked by a digital filter.
The advantages of performing a digital filter simulation include the
following. First, researchers can obtain data without performing expensive and
time-consuming experiments. Second, the simulation provides a flexible means
of varying the input and system parameters and observing the effect of these
changes. A third advantage of simulations is that any desired signal in the
model can be observed. In the physical system, it may be very costly or even
impossible to observe the value of certain parameters that are important in the
operation of the system. The following cases illustrate these advantages.
Simulating a computer communication system. One common digital data
communication system that allows computers to transfer data uses the
frequency-shift keying (FSK) scheme. In the FSK system, the binary digits, 1
or 0, are transformed into sections of a sinusoidal waveform. As shown in Fig.
1.10, the frequency of a particular section is determined by the binary value to
be transmitted. A simulation program including a digital filter to model the
transmitter, detector and the transmission channel may attempt to answer the
following questions: What frequencies should be used? How long should the
duration of each segment be? What detector is suitably reliable but reasonably
expensive or fast? To investigate the complicated effects of interference,
random noise can be added to the transmitted signal and the resulting
characteristics of a given detector can be examined. All these procedures can
be performed before any physical equipment is manufactured or installed.
Such a simulation program is developed by the student in the projects at
the end of the chapters. In one project, the transmitter is simulated by a
program that generates the waveform that corresponds to a given binary
sequence. Other projects consider different digital filters that operate on the
received sinusoidal sequence to retrieve the value of the binary digit
transmitted.
Simulating a robot arm. Another advantage of simulation programs is that the
various parameters acting in a physical system can each be examined in detail.
For example, consider the simulation of a robot arm, shown in Fig. 1.11. The
actual arm is constructed of some material that has mass and elasticity. The
input to such an arm is shown as a rotational displacement that is modeled as a
step function at one joint of the arm. The desired output may be the location
of the manipulator at the end of the arm. The time response of the output may
be important for ensuring that the arm construction satisfies the specification
for quickness or overshoot.
One simple digital filter model of the arm is a second-order filter whose
coefficient values are related to the mass and elasticity of the material. The
physical input signal is converted into an input sequence, as shown in the
figure, and the output sequence of the filter indicates the response to be

INTRODUCTION

(a)

11

--"0"---,

r--"O"

"t"----i>I1IE-E

x(t)

~----~-----+--~--~-+--~----~-----,time

(b)

ftGURE 1-10
Frequency-shift keying (FSK) signals. a) Example of continuous-time signal from actual system. b)
Simuated sequence produced by computer.

{e(D)} = {

8(D)~

{e(n)}
___ ...r_.:::_..... Ili n)}
I

-----~:--;..-

= (y(n)}

:111 III II III I 1':,

o

{lin)}

~tlill [ rill I "":,

o

!!

ftGURE 1-11
Simulating a robot a'lll by applying an excitation sequence and computmg the response sequenct'

INTRODUcnON TO DIGITAL SIGNAL PROCESSING

expected if the model is accurate. The value of this simulation model is that
the coefficient values of the digital filter can be adjusted to produce the desired
response. Then the final values can be employed to determine the physical
dimensions of the robot arm.

1.S REPLACING ANALOG FILTERS

WITII DIGITAL FILTERS
One of the first applications of digital filters was to replace analog filters,
composed of resistors R, capacitors C and inductors L. This replacement is
usually performed to overcome some of the inherent limitations of analog
components, including the fluctuation of the component values with temperature and age. Another disadvantage of analog components is their large
physical size, especially in the case of large-valued inductors and capacitors.
A further disadvantage of analog filters is the difficulty encountered in
changing the component values iq an existing filter. In many cases, the
equivalent digital filter is implemented as a sequence of instructions in a
computer, which can easily be changed with an editor, while soldering may be
required to change a physical component in the analog filter. When the digital
filter is implemented in hardware, the multiplier coefficients are often stored in
a read-only memory (ROM), which can readily be replaced if a filter
modification is desired. These coefficients can also be stored in read/write
memory and be changed automatically as the nature of the input signal
sequence changes. This added flexibility of digital filters has created an area of
research in so-called adaptive filters that is currently very active.
Replacing an analog filter with a digital filter in a given application
involves several considerations. For illustration, let us examine the situation in
Fig. 1.12, which shows a continuous-time input signal x(t), a simple analog
RC lowpass filter and the output signal y(t). To replace this analog filter with
the equivalent digital filter, analog-to-digital conversion (ADC) must be
performed to produce a discrete-time input sequence {x(n)} from x(t). The
input sequence {x(n)} is then processed with a digital filter, implemented with
special-purpose hardware or as a computer program, to produce the output
sequence {y(n)}. Finally, a digital-to-analog conversion (DAC) must be
performed to reconstruct the continuous-time signal y(t). For simple analog
filters, the added expense involved in the ADC, DAC and the hardware to
implement the digital filter does not usually make this substitution economically feasible. But with the proliferation of microcomputers and the reduced
cost of analog-to-digital converters, performing a filter operation with digital
circuits is becoming common.
In some applications, the analog filters themselves may be ~tly. In
geophysical and biomedical applications, signals with very low-frequency
components, typically less that one Hertz, require very bulky and expensive
inductors and capacitors. In these cases, digital filters are commonly employed
as replacements. In some applications, such as those involving very high

INTRODUCTION

13

ANALOG FILTER
x(t)

-"'I)

DIGITAL FILTER

ftGURE lU
Replacing an analog filter with a digital filter requires the additional operations of an
analog-ta-digital conversion (ADC) and a digital-to-analog conversion (DAC).

frequencies, analog filters still have the advantage. But with decreasing cost
and increasing speed of digital circuits, more and more analog filters are being
replaced with their digital counterparts.

1.6 OVERVIEW OF THE COURSE

The digital filter design and signal processing concepts described in this book
can be divided into five categories: analysis, synthesis, the fast Fourier
transform, finite-accuracy effects, and inverse filters.
ANALYSIS. We analyze digital signal processing algorithms by determining
their time and frequency domain characteristics. This allows us to attain an
intuitive understanding of their capabilities and limitations. The time-domain
performance of digital filters is described in terms of the filter's unit-sample

14

INTRODUCTION TO DIGITAL SIGNAL PROCESSING

response sequence, denoted by {h(n)}. This sequence is analogous to the

impulse response of analog filters. A convolutional equation allows us to
determine the output sequence {y(n)} from the input sequence {x(n)} and the
unit-sample response {h(n)}. The unit-sample response sequence also permits
us to determine whether a filter is stable. Linear difference equations provide
an alternate time-domain description that is be useful for implementing digital
filter structures.
Filter specifications are commonly dpressed in the frequency domain.
Hence, we examine the relevant aspects of the Fourier domain that define the
frequency domain properties of signals and filters. The Fourier transform of
the unit-sample response {h(n)} is the transfer function H(e iw ) of the filter and
it describes the gain of the filter at different frequencies. The Fourier transform
of a data sequence {x(n)} is called the spectrum X(e iw ) and it defines the
frequency content of the signal. We explore the problems encountered when
generating a discrete-time sequence by sampling a continuous-time signal and
when reconstructing the continuous-time signal from the discrete-time sequence. The discrete Fourier transform (DFT) allows us to evaluate the
Fourier transform by computer. The consequences of employing' the DFT for
signal processing and spectral analysis applications are investigated.
An alternate and more general analytic technique than the Fourier
transform is the z-transform. The system function H(z) is defined as the
z-transform of the unit-sample response {h(n)} and is used for the analysis and
synthesis of digital filters. The complex z plane is the domain of interest for the
z-transform. The representation of a digital filter as a collection of poles and
zeros in the z plane will provide a useful interpretation of the filter frequency
response.
SYNTHESIS. After having a firm and intuitive understanding of the time and

frequency behavior of digital filters, the goal is then to design and implement a
digital filter to perform a desired task. The synthesis procedure involves
combining the digital filter elements and specifying the values of the multiplier
coefficients to satisfy a desired specification.
Most often, the transformation that is to be performed by the digital filter
is specified in terms of the magnitude response IH(eiw)l. The magnitude
response for an ideallowpass filter, for which the frequency components in the
input sequence from zero to some cutoff frequency Wc are passed through the
filter unaffected, while all th~ others are completely eliminated from the output
sequence, is shown in Fig. 1.13(a). This ideal cannot be achieved in practice,
but must be approximated. A practical specification that can be satisfied by a
digital filter is shown in Fig. 1.13(b). This desired specification allows slight
deviations from the ideal response: E in the passband, () in the stopband, and a
nonzero frequency interval for the transition between the passband and
stopband. A successful digital filter design produces a magnitude response that
falls within the allowed regions.
We will consider two different digital filter types to satisfy a given

INTRODUCTION

I<

PASS BAND

0
(b)

=-jE

STOPBAND

w,

;.

15

IHdc>"cd(eJ"')1

Successful design

1-

I
I
I
I
I
I

-----------1-----------~.~~--~~--~
OL-------~--------~~~=--=~~

I<

PASS BAND

~IE
wp

TRANSITION
RANGE

STOP BAND

~I W

Ws

FlGURE 1-13
Magnitude response curves for a digital low pass filter. (a) Ideal low pass rlter response. (b)
Practical low pass filter specification allows derivations from the ideal response. A successful
digital filter has a magnitude response that falls within the allowed regions.

specification: infinite impulse response (IIR) filters and finite impulse response
(FIR) filters. IIR digital filters are usually designed by extending classical
analog filter design procedures. An I1R filter design is typically accomplished
in three steps. First, an analog lowpass filter is designed to meet the desired
passband specification. The most commonly employed are the Butterworth,
Chebyshev and elliptic analog lowpass filter design procedures. Second, an
analog-to-digital filter transformation is employed to obtain a digital lowpass
filter. The impulse-invariance method and the bilinear z-transform method are
described for transferring these analog designs into their digital counterparts.
In the final step, a frequency transformation is employed to obtain a highpass,
bandpass or bandreject filter.
/
The design procedures for FIR filters have been develope~SpeciallY for

16

INTRODUCTION TO DIGITAL SIGNAL PROCESSING

digital filters and do not have analog filter counterparts. Two FIR filter
structures will be investigated. The first involves implementing the desired
filter unit-sample response in a tapped-delay line structure. Since the desired
response often has an infinite number of nonzero elements, a finite-duration
window must be applied to eliminate all but a finite number of these elements.
We investigate the effect of this truncation on the resulting filter's transfer
function. The second structure is the frequency-sampling structure, which
employs a comb-filter followed by a bank of resonators and implements the
filter by specifying the magnitude response at a finite number of frequency
points. A comparison of the two design procedures is performed to determine
under which conditions one is more efficient than the other.
FAST FOURIER TRANSFORM (FFf). In addition to digital filters, we also
explore the use of the fast Fourier transform (FFT) algorithm for processing
data and performing spectral analysis. The FFT is an elegant and efficient
method for computing the discrete Fourier transform that allows the filtering
operation to be performed in the frequency domain. We present both an
intuitive description of the FFT, as well as a more rigorous mathematical one,
to indicate how this efficiency is achieved. The utility of employing the FFT for
filtering applications is described and compared with the methods of direct
digital filter implementation.
FINITE-ACCURACY EFFECTS. An important concern for implementing practical digital filters involves finite-accuracy effects caused by the limited
precision to which the sequence and coefficients values are represented within
the computer. Similar errors are also introduced in the numerical round-off
routinely performed in tedious hand calculations. In the analysis and synthesis
discussions above, we dealt with the ideal, or infinite-precision, case. To
process data on a computer, the data values must first be quantized to a finite
number of bits, usually equal to or smaller than the word size on the computer.
The effects of this quantization error will be investigated. The coefficients used
in the digital filter multipliers must also be quantized. The slight change in
these coefficient values modifies the frequency response of the filter. A third
source of error occurs when several numbers are added or multiplied and the
result excee'ds the limits of the computer. In this case, the result must be scaled
and truncated, introducing still another source of error that can produce
unwanted oscillations at the output of some IIR filters. We consider analytic
methods of characterizing these three common quantization problems.
INVERSE FILTERS. In the final chapter, we discuss the interesting and
important topic in signal processing called inverse filtering. This topic integrates many of the digital signal processing concepts presented in the previous
chapters to solve the problem of removing distortions from observed signals
that are introduced by measuring instruments or by physical transmission
media. It is hoped that this topic whets the student's appetite for pursuing
further study in this most interesting area of research and development.

JNTRODUcnON

17

1.7 DESCRIPTION OF THE COMPUfER

PROJECfS
Along with the standard problem sets at the end of each chapter, a set of
projects to be performed on a computer is also included. The pedagogical
value of these projects cannot be overstressed, since they are indicative of the
problems encountered in practice. The computer is a powerful tool for signal
processing and the projects given in this book illustrate its flexibility. Most of
the projects in Chapters 2 through 6 can be analytically determined with little
effort. The analytic results are meant to verify the computer results and to
instill some faith in the algorithms. As the analysis becomes more difficult, the
computer results can help guide the analytic procedure. The design projects in
Chapters 7 through 11 are almost impossible to do analytically, as are most
practical filter designs. By that stage in the course, the intellectual insights are
achieved almost exclusively from the computer results. Using the computer in
this fashion is the basis for computer simulations, which are often performed to
explore the performance of complicated systems before they are implemented.
The subroutines for the projects are meant to be written by the st\1dent,
in FORTRAN, the current universal language for signal processing, or in the
increasingly popular C or PASCAL (BASIC is slightly less desirable, because
of its slow speed, but acceptable). To minimize the amount of programming, a
structured approach is used in which the routines written for the initial projects
are used again in subsequent projects. Each routine performs a simple, and
hence easily verifiable, operation. Constructing more sophisticated programs
from these simpler building blocks. aids the student in conceptualizing the
digital filter operation a!,ld for trouble-shooting the more complicated projects
later in the course. By the end of the course, the student should have an.
intuitive understanding of digital signal processing that is invaluable for
applying it to novel situations. As an added benefit, students will also have a
flexible library of signal processing programs that they have written themselves, and hence understand intuitively.
Graphic outputs produced by the projects are displayed on the terminal
and on a printer, being generated by programs provided in Appendix A. These
plots, although having crude resolution, are adequate for illustrating the ideas.
They also have the advantage of not requiring special plotting hardware. The
plotting subprograms can also easily be modified to display the results on other
types of digital plotters and display terminals.
The experience gained through this course should whet the appetite for
the study of advanced topics in filter design and signal processing. References
to material describing advanced topics are provided at the end of each chapter.
The student will also be well prepared to apply digital filtering techniques to
practical problems.
After this course, the student can easily become facile with the various
programs that are available for automatically designing optimal digital filters.
Examples are the digital filter program packages offered by the IEEE, or in
SIG, the software package from the Lawrence Livermore National Laboratory

18

lNTRODUcnON TO DIGITAL SIGNAL PROCESSING

at the University of California. An important prerequisite for employing these

programs successfully is the intuitive understanding of the tradeoffs and
subtlety involved in the methods for digital signal processing. This is the
material stressed in this book.

1.8 SUMMARY
In this chapter, we have provided the motivation for digital filtering, defined
the components of a digital filter, and indicated some common applications.
An overview of the course material has described the inter-relationship of the
topics that describe the theory and application of digital filters to signal
processing.

FURmER READING
IEEE Acoustics, Speech and Signal Processing Magazine. This bimonthly magazine of the
Acoustics, Speech and Signal Processing Society of the Institute of Electronics and Electrical
Engineers (IEEE) provides tutorial articles on various topics covering digital signal processing
and digital filters.
IEEE Transactions on Acoustics, Speech and Signal Prl?cessing. This bimonthly journal publishes
articles describing the most recent results in the development and implementation of digital
filters and digital signal processing techniques.
Rabiner, L..R. and C. M. Rader (eds.): Digital Signal Processing, IEEE Press, 1972; and Digital
Signal Processing ll, IEEE Press, 1975. [These books contain reprints from various journals
that provide an organized and historical treatment of digital filter theory and applications.]

REFERENCES TO TOPICS
Pascal
Cooper, Doug, and Michael Clancey: Oh! Pascal!, W. W. Norton, 1982.
Jensen, Kathleen, and Niklaus Wirth: Pascal User Manual and Report, Springer-Verlag, Berlin,
1985.
Mallozzi, John, and Nicholas De Lillo: Computability with Pascal, Prentice-Hall, Englewood
Cliffs, NJ, 1984.

C language
Harbison, Samuel, and Guy Steele Jr.: C, A Reference Manual, Prentice-Hall, Englewood Cliffs,
NJ,1984.

FORTRAN
Didday, Richard, and Rex Page: FORTRAN for Humans, West, 1981.
McCracken, Daniel: FORTRAN, Wiley, New York, 1967.
BASIC

Maratek, Samual: BASIC, Academic Press, New York, 1975.

INTROOycnON

19

Comptder (OIIUIIlUIicatioas
Schwartz, Mischa: Telecommunication Networks, Addison-Wesley, Reading, Mass. 1987.
Tanenbaum, Andrew: Computer Networks, Prentice-Hall, Englewood Cliffs, NJ, 1981.
Lathi, B. P.: Modem Digital and Analog Communication Systems, Holt, New York, 1983.

DIPDI sIpaI processIq propam packages

Programs for Digital Signal Processing, IEEE Press, New York, 1979.
Lager, Darrel L., and Stephen G. Azevedo: SIG -A General-Purpose Signal Processing Program.
Proceedings of the IEEE, Vol. 75, No.9, pp. 1322-1332, 1987.

RoIIot

_aIadoa

IEEE Transactions on Pattern Analysis and Machine Intelligence.

Speedt proceaiq

IEEE Transactions on Acoustics, Speech and Signal Processing.

MecIkaI sIpaI procesIIbIc

IEEE Transactions on Biomedicfll Engineering.

CHAPTER

2
DISCRETE-TIME
DESCRIPTION
OF SIGNALS
AND SYSTEMS

,r

2.1 INfRODUcnON

. In this chapter we investigate procedures to describe the time-domain behavior

of discrete-time systems, of which digital filters are a special class. After
establishing the discrete-time signal notation, we define several sequences that
are important for the llnalysis of digital filters. The superposition principle for
linear systems is then employed to derive useful relationships governing digital
filters. The unit-sample response is presented as a time-domain descriptor of a
discrete-time system. The linear convolution operation relates the output
sequence of the filter to the input sequence and the unit-sample response.
Linear constant-coefficient difference equations are employed for the analysis
and synthesis of digital filters. A procedure to write a digital filter program
directly from these time-domain equations is described.

2.2

DISCRETE~TIME

SEQUENCES

In this book, we will be dealing with discrete-time sequences, which will be

denoted either by {x(n)}, or by x(n); for -Nl $;n $;N2, where Nl and N2 may
possibly be infinite. As mentioned in Chapter 1, the term discrete-time can also
refer to sequences of other physical parameters, such as discrete-location for

DISCRETE-TIME DESCRIPTION OF SIGNALS AND SYSTEMS

21

image processing applications, discrete-distance for sonar ranging systems, or

discrete-position for robotic arm sensors_ The particular value of the sequence
at time k will be denoted simply by x(k), with no additional qualifiers on
the index. In this book, we are primarily concerned with sequences that are
produced by sampling a function of a continuous-valued variable at equally
spaced intervals. These, in fact, are the most common types of sequences in
digital filter applications. We now consider several important sequences that
are used throughout this book for analyzing digital filters.
UNIT-8AMPLE SEQUENCE. The unit-sample sequence contains only one
nonzero valued element and is defined by

d n = { 1 for n = 0

()

0 otherwise

(2.1)

This sequence is shown in Fig. 2.1 and will play an important role as the input
sequence to a digital filter. Since the unit-sample sequence contains only one
nonzero element, the form of the output sequence then directly indicates the
effect produced by the digital filter.
The delayed unit-sample sequence, denoted by {d(n - k)}, has its
nonzero element at sample time k:
d(n - k)

= {~

for n = k
otherwise

(2.2)

The unit-sample sequence is also an important analytic tool because of its

sifting property, by which a particular element can be extracted from the entire
sequence. To demonstrate this property, we note that the value of {x(n)} at
time k can be written as
DO

x(k)= ~ d(n-k)x(n)

(2.3)

n=-OCI

as shown in Fig. 2.2. This follows because d(n - k) is nonzero only when its
argument is zero, or when n = k. When multiplying the two sequences together
term by term, only the single element x(k) remains nonzero and the sum
evaluates to x(k). Hence, the unit-sample sequence can be used to select, or
sift, a single. element from the entire sequence. This sifting property will be
important below in the derivation of the convolution relationship between the
digital filter input and output sequences. A more convenient, but equivalent,
form of Eq. (2.3) that is used below employs the time-reversed version of the
{d(n)}

... r. ..

-2-10 123 ...

FIGURE 2-1
Unit-sample sequence.

22

INTRODUCTION TO DIGITAL SIGNAL PROCESSING

(x(n)}

... 111111!]11
I
...

k ... -I 0 1 2

:>

(d(n-k)}

.c.
k

:>

n
(u(n)}

... ...IIIIIII
o

FIGURE 2-2
Sifting property of the unit-sample sequence.

FIGURE 2-3
Unit-step sequence.

unit-sample sequence having its nonzero element at n

= k:

00

x(k)=

d(k-n)x(n)

(2.4)

n=-CIO

UNITSTEP SEQUENCE. The unit-step sequence is defined as

u(n) = {~

forn~O

otherwise

(2_5)

and is shown in Fig. 2.3. The unit-step sequence is used for defining the
starting point of a sequence in analytic expressions. For example, the sequence
an
for n ~O
(2.6)
x(n) = {
otherwise

0,'

can be written simply as

for aU n

(2.7)

SINUSOIDAL SEQUENCES. Sinusoidal sequences play an important role in

the frequency-domain analysis of digital filters. Given some radian frequency

w, we can form the discrete-time sinusoidal sequences {cos(wn)} and
{sin(wn)} by evaluating these functions at the corresponding value of the
argument. Examples are shown in Fig. 2.4. Strictly speaking, the units of ware
radians per sampling interoal. By convention, the sampling interval is taken as
a dimensionless quantity, allowing OJ to be specified in radians. In the next
chapter, we will consider the sampling of continuous-time signals and relate OJ
to continuous-time frequency values.

DISCRETETIME DESCRImON OF SIGNALS AND SYSTEMS

(cos(wn)}

"'1.

1II

w=.~

. iiI.
-I 0 1

(sin(wn)}

1I

23

....
21 II

(o)=~

...
tIllllt
10

III1

-2- l o

1 2

FIGURE l-4
Two discrete-time sinusoidal sequences.

Discrete-time sinusoids exhibit an interesting property that is important

for evaluating the Fourier transform of a discrete-time sequence: The set of all
distinct values that a discrete-time sinusoidal sequence can have occurs for
values of co lying in the range -:re:s; co <:re.
To illustrate this property, let us consider two frequency values, Wo is in
the interval [-:re, :re), and WI outside this range. For example, let
COl

= Wo + 2m:re

(2.8)

where m is a nonzero integer. The values of the cosine sequence at frequency

COl are given by
{cos(co1n)} = {coscoo + 2m:re)n)}

= {cos(won)}

(2.9)

This last result follows because the cosine function is periodic with period 2:re
and the argument is won plus an integral number ( =tnn) of 2:re. In Fig. 2.5,
this effect is shown for COo = 0.1:tr and COl = 2.1:re. This result is also true for the
sine sequence.
COMPLEX EXPONENTIAL SEQUENCE. Before we consider the complex
exponential sequence, we will establish the notation for complex numbers and
the complex plane used throughout this book. Let c be a complex number.
Then it can be written as
c =a + jb

(2.10)

where a is the real part, b is the imaginary part, and j = (_1)112. The
representation of c on the complex plane is shown in Fig. 2.6. An alternate
form of a complex number is the complex vector notation:

c = lei ei Arg[c)
where

(2.11)

lei is the magnitude and Arg[c] is the argument or phase, which is

24

INTRODUcnON TO DIGITAL SIGNAL PROCESSING

{cos (0. lIT) }

...
{cos (2.IIT)}

FIGURE 25
Two analog signals having different frequencies can proouce the same values for a discrete-time
sinusoidal seque\lce.

defined by the angle with respect to the real axis in the complex plane. The
following geometric relationships can be derived from Fig. 2.6:
and

Icl 2 = a2 + b2

(2.12)

Arg[c] = arctan[bla].

(2.13)

The complex exponential sequence is important for a compact frequency

analysis of digital filters, and is given by
{e("+i",)n} = {e<7n ei"'n}
(2.14)
When both a and ro are real numbers, each element of the sequence is
expressed in complex vector notation with magnitude e"" and phase ron.
A special case of this sequence occurs when~ a = O. The exponential
sequence can then be related to the sinusoidal sequences by the geometric
Imaginary
c=a+jb

-------~---'L-.----I:--~

Real

FIGURE 2-6
The complex plane.

DISCRETE.TIME DESCRIPTION OF SIGNALS AND SYSTEMS

2S

Imaginary
c=elNf'I

------4<--"---'--_ Real

Imaginary

-----~-+---_:_.."..__:_~

Real

Imaginary
2)sin(wn)

-e-"'"

eJ"'"

____--l..-______

Real

F1GURE 27
Interpretation of the Euler identities in the
complex plane.

relationships shown in Fig. 2.7, which are known as the Euler identities. These
are given by
ei",n = cos( ron) + j sin( ron)
(2.15)
ei",n + e-i",n
(2.16)
cos(ron) =
2
sin(ron) =

ei",n - e-i",n
2j

(2.17)

The property of sinusoidal sequences, that all values can be represented

with ro in the range [-Jr, Jr), is also true for complex exponentials. To
illustrate this, consider rol = roo + 2mJr, where m is any integer. Then
(2.18)

26

1NTR0DUcnON TO DIGITAL SIGNAL PROCESSING

since ej2mn = 1, when m is any integer. This identity can also be demonstrated
in the complex plane. Increasing the phase of a complex number by a is
equivalent to rotating the complex vector by angle a, counter-clockwise if
a> 0 and clockwise if a < O. If a = 2m:r, where m is any integer, then the
complex vector makes m complete revolutions about the origin and returns to
its original orientation.

2.3 SUPERPOSmON PRINCIPLE FOR

LINEAR SYSTEMS
The digital filters considered in this book are a special class in the general
category of linear systems. Any discrete-time system can be thought of as
transforming an input sequence {x(n)} into an output sequence {y(n)}
through a transformation denoted by G { . }. However, two characteristics
define a linear system: (1) for a given input sequence, a change in the
amplitude scale of the input sequence results in the same amplitude scale
change in the output sequence, and (2) if two input sequences are added and
the sum is applied to the system, the resulting output is the same as the sum of
the responses to the individual inputs. These ideas can be formalized in the
statement of the superposition principle.
SUPERPOSmON PRINCIPLE. Given a system producing a transformation
G{ .}, two arbitrary input sequences {xI(n)} and {x2(n)} that result in
{YI(n)} = G{xI(n)} and {Y2(n)} = G{x2(n)}, and two constants a and (3; then

the system is defined to be linear when

G{axl(n) + (3x2(n)}

= aG{xI(n)} + (3G{x2(n)}
= a{YI(n)}

+ (3{Y2(n)}

(2.19)

The following example illustrates the applicatiOn of the superposition

principle.
Example 2.1. Linearity of the 3-sample averager. The transformation produced
by the 3-sample averager is given by

y(n) = Ux(n + 1) +x(n) +x(n -1)] = G{x(n)}

To test the linearity, we apply the scaled sum of two arbitrary sequences, or
{axt(n) + (3x2(n)}, to the averager. The output is then equal to

G{axt(n) + (3x2(n)} = Uaxt(n + 1) + (3x2(n + 1) + axt(n) + (3x2(n)

+ axt(n -1) + f3x2(n -1)]

="3a [xt(n + 1) + xt(n) + xt(n -1)]

+~[X2(n + 1) +x2(n) +x2(n -1)]

= a{Yt(n)}

+ (3{y2(n)}

DISCRETETlME DESCRIPTION OF SIGNALS AND SYSTEMS

27

This last equation is simply the sum of the two outputs that would have been
observed if only {crxt(n)} or {,8x2(n)} had been applied. Hence, the 3-sample
averager is a linear system.

Example 2.2. A nonlinear system. Let us consider a square-law device, in

which the output is the square of the input:
{y(n)} = G{x(n)} = {x 2(n)}
By applying two scaled inputs into this system, we have at the output:
G{crxt(n) + ,8x2(n)} = [crxt(n) + ,8x2(nW
= ~xi(n)

+ ,82x~(n) + 2a,8Xt(n)x2(n)

which is not equal to ~xi(n) + ,82x~(n), the output required for the system to be
linear. Hence, the square-law device is a nonlinear system.

2.4 UNIT-SAMPLE RESPONSE

SEQUENCE
One important consequence of the superposition principle is that it allows us to
derive a convenient relationship between the input and output of a linear
system in the time domain. Applying the sifting property-of the unit-sample
sequence given by Eq. (2.4), we can write
{y(n)}

= G{x(n)} = GL~oo d(n = G{ ..

k)X(k)}

+ x( -1)d(n + 1) + x(O)d(n) +x(l)d(n -1) + ... } (2.20)

Since the system operates through its transformation G{ . } on sequences,

we can view the output of a system as the superposition of responses to
delayed and scaled unit-sample sequences, as shown in Fig. 2.8. The
unit-sample sequence having its nonzero element at n = k, d(n - k), is scaled
by x(k). Applying the superposition principle, we can 'write the output as
{y(n)}

= ... + x( -1) G{d(n + I)} + x(O) G{d(n)} + x(l) G{d(n -I)} + ...
00

= 2:

x(k) G{d(n - k)}

(2.21)

k=-oo

G{d(n - k)} is the output sequence produced by the system when the input is
a unit-sample sequence having fts nonzero element at n = k. This output is
called the unit-sample response, and is sufficiently important to merit its own
symbol:
G{d(n - k)}

= {h(n, k)}.

(2.22)

In words, {h(n, k)} is the response of a linear system whose input is a

unit-sample sequence with nonzero element at position n = k. The following
examples illustrate the concept of the unit-sample response.

28

INTRODUcnON TO DIGITAL SIGNAL PROCESSING

(..)}
If:
{.t(II)}

---1

~(y(.)}

LINEAR
SYSTEM

= {d(lI)}

I..
0

{Y(II)}

>
~

II

= {h(lI)}

.. IIIt ...
-I 0 I

II

Then:
{Y(II)}

{.t(II)}

~I

0 I

II

-I

.t(

)(d(lI+ I)}

-I
L~-I:

II

o
---.........l-Tl"""'
r-r......'""'...........;..~I .t(-l){h(lI+ I)}
+

.t(O){d(lI)}

+
.t(I){d(lI+ I)}

I.t!O~
0

II

T , >
o

.t(O){h(lI)}

II

... r.
oI

II

Ji'lGURE U
The input sequence to a \inear system can be viewed as the superposition in time of a set of scaled
unit-sample sequences. The output sequence is then the superposition of the scaled unit-sample
responses.

Eumple 2.3. Uait-sample response of the 3-sampIe .veneer. Let

y(n) =lx(n + 1) + lx(n) + lx(n -1)

To find the unit-sample response {h(n, O}}, we set {x(n}} = {d(n}} and solve for
the values of y(n) for each n:
For n = -2,
y( -2) = id( -1) + Id( -2) + ld( -3) = O.
Forn S -2,
y(n) = 0, since d(n) = 0 for n <0.

DISCRETETlME DESCRIPTION OF SIGNALS AND SYSI'EMS

29

{h(n,O) }

.:111 ..

3 21 0 1 2 3

Forn=-I,
For n =0,
For n = 1,
For n=2,

FIGURE 2-9
Unitsample response
sample averager.

y(-I) = id(O) + !d(-I) + id(-2) =

of nonrecursive three-

i.

y(O) = !d(I) + i~(O) + id( -1) '" }.

y(I) = id(2) + id(I) + id(O) = ...
y(2) = id(3) + id(2) + id l" . O.

For n>2, y(n) =0, since the nonzero value of {d(n)} has moved out of the
memory of this filter. Since {d(n)} is the input, the output is defined to be
{h(n, O)}. Hence
h(n, 0) =

{~,

for-I""n",,1
otherwlse

This sequence is shown in Fig. 2.9.

The unit-sample response is computed by starting at a point in time, n,

for which all the previous inputs and the contents of the delays and advances
are all equal to zero. This state of the filter is known as the zero initial
condition.
EumpIe 2.4.

Unit-5lUllple

respoIIIe

of. first-order recursive filter. Let

y(n) = {~(n -1) +x(n)

for n 2:0
otherwise.

To find {h(n, O)}, we let {x(n)} = {d(n)} and apply the zero initial condition.
For n<O, y(n) =0, because d(n) is zero for n <0 and y(-I) =0.
For n = 0, y(O) = ay( -1) + d(O) = 1.
For n = 1, y(I) = ay(O) + d(I) = a.
For n = 2, y(2) = ay(I) + d(2) = a2
For n =k,.y(k) =ay(k -1) + d(k) = a".
Hence,
h(n, 0) =a"u(,n)

for all n

Some typical sequences are shown for different values of a in Fig. 2.10.

This concept of a unit-sample response can also be applied to timevarying systems, in which the multiplier coefficients vary with n, as illustrated
in the next example.

30

INTRODUCfION TO DIGITAL SIGNAL PROCESSING

h(n,O)

II

for all n

.'.
0<.<1

..1111 J
o

= a"u(n)

t , " ";"

)0

-I <a<O

a<-l

FIGURE 2-10
Four possible unit-sample response sequences that can be produced by first-order recursive
filter.

Example 2.5. Unit-sample response of a time-varying system. Let

Y()
n =

{ o~IY(n-l)+x(n)
n+

forn""O

otherwise.

Applying the zero initial condition and using the same approach as in the
previous example, we have
h(O, 0) = (1)(0) + 1 = 1

h(I, 0) = 0>(1) + 0 = ~
h(2, 0) =
h(3, 0)

<no> = ~

= mm = -i4

If the unit-sample is applied at time n = 1, i.e. {x(nn = {d(n

{y(n}} = {h(n, In, where
h(O, 1) = (1)(0) + 0 = 0

-In, then

h(l, 1) = G)(O) + 1 = 1

h(2,I)=m(I)+O=i
h(3,I)=mm=fi

In

On.

Note that {h(n,

is a different-valued sequence from {h(n,
This example
illustrates that, for a time-varying filter (whose coefficient values vary with time
n), the values of the unit-sample response depend upon the time of application of
the unit-sample sequence.

DISCRETE-TIME DESCRIPTION OF SIGNALS AND SYSTEMS

31

The general time-domain relationship between the system output value at

time n and the input values can be obtained by combing Eqs_ (2-21) and
(2.22), to give

y(n) =

(2.23)

x(k) h(n, k)

k=-oo

This form is called a convoLutionaL reLationship. This equation can be simplified

when it is applied to time-invariant systems, which are discussed in the next
section. The references describe adaptive filters whose coefficient values
change to meet various conditions.

2.5 TIME-INVARIANT SYSTEMS

We will further restrict our attention in this book to the important class of
Linear time-invariant (LTI) systems, or those whose transformations do not
change with time. In simplest terms, for an LTI system, a shift of no samples in
the input produces onLy a corresponding shift in the output. That is, if
{yen)} = G{x(n)}, then

{yen - no)}

= G{x(n -

no)}

(2.24)

for all no

One characteristic of an LTI system is that the form of the unit-sample

response does not change with the time of application of the unit-sample
sequence. Then, rather than being a function of two variables, nand k, the
unit-sample response becomes a function of only n - k, or the time as
measured from the application of the unit-sample. Analytically, we then write

G{d(n - k)}

= {hen, k)} = {hen -

(2.25)

k)}

This equation implies that a shift in the time that the unit-sample sequence is
applied to the system results only in a corresponding time shift in its
unit-sample response, as shown in Fig. 2.11. A time-invariant digital filter can
easily be recognized, since its coefficient values do not vary with time.
{d(n))

.'I

0 1 2

.. 1.

o 1

.. IJ
o1

Tt
2

(h(n.k)} = {h(n-k)}

(h(n))

, , ";'
n

.... IITt, ..;

o

F1GURE 211
Shift-invariance properly of the unit-sample response for a time-invariant system.

32

INTRODUCTION TO DIGITAL SIGNAL PROCESSING

Applying Eq. (2.25), the convolution relationship of a LTI system can be

written as
00

y(n)=

x(k)h(n-k)

(2.26)

= {x(n)} * {hen)}

(2.27)

k=-cc

or, more concisely,

{yen)}

where * will denote the convolution operation. In evaluating a convolution

expression, it is often helpful to note that convolution is a commutative
operation, i.e.
{yen)} = {x(n)}

To demonstrate this, we let m

yen) =

L
m=cc

* {hen)} = {hen)} * {x(n)}

=n -

(2.28)

kin Eq. (2.26) to get

00

x(n - m) hem)

= L

h(m)x(n - m)

(2.29)

m=-oo

The desired result is obtained when m is replaced by k. The limits of the first
sum in Eq. (2.29) appear in reverse order, but can be reordered without
consequence because the same terms appear in the evaluation of the sum.
The convolution operation is also distributive, i.e.
{hen)}

* {xt(n) +x2(n)} = {hen)} * {xt(n)} + {hen)} * {x2(n)}

(2.30)

This result can be shown directly by applying the superposition principle.

COMPUTING A DISCRETE-TIME CONVOLUTION. Since the convolution

operation is one of the most common and important in digital signal

processing, we present a step-by-step procedure for the evaluation of (2.26).
Step 1. Choose an initial value of n, the starting time for evaluating the
output sequence. If {x(n)} starts at n = nx and {hen)} starts at n = nh, then
n = nx + nh is a good choice. The values of x(n) for n < nx and of h(n) for
n < nh are then assumed to be zero.
Step 2. Both sequences are expressed in terms of the index k. Note that
x(k) is identical to x(n), while hen - k) is a time-reversed version of hen). If
hen) starts at n = nh' thea the last point of hen - k) occurs at k = n - nh.
Step 3. The two sequences are multiplied element by element and the
products are accumulated over all values of k. This sum of products then
becomes the value of the output for index n, or y(n).
Step 4. The index n is incremented and steps 2, 3 and 4 are repeated until
the sum of products computed in Step 3 is zero for all the remaining values
of n.

DISCRETE-TIME DESCRIPTION OF SIGNALS AND SYSTEMS

33

These four steps are illustrated in the following examples.

Example 2.6. Convolution of two finite-duration sequences. Let

x(n) =

{~

for -1 sn s 1
otherwise

h(n) =

{~

for -1 Sn S 1
otherwise

and

The evaluation of {y(n)} = {h(n)}

given above is shown in Fig. 2.12.

* {x(n)}

by Eq. (2.26) using the four steps

Step 1. The evaluation will start at n = -2.

Step 2 Both sequences are shown indexed by k, with {h(n - k)} being the
time-reversed version of {h(n)}. Since {h(n)} starts at n = -1, the last point of
{h(n -k)} occurs at k =n + 1. When n = -2, k =-1.
Step 3. Computing the product of the pair of sequence values for each k,
only the k = -1 terms yield a nonzero product. This result then becomes the
value of y(-2).
Step 4. The value of n is incremented and the process above is repeated
until n = 3. The relative p'ositions of the two sequences are shown for the values
of n that produce nonzero values for y(n).
Example 2.7. Convolution of a 6nite-duration seqnence with an infiniteduration sequence. Let
forOSn s2
otherwise

x(n)={n+1,
0,
and

h(n) = a"u(n)

for all n

We will perform the convolutio.n in two ways. First, using Eq. (2.29), we have
~

y(n)=

2:

2:

aku(k)x(n-k)

for n -k =0,
for n -k = 1,
for n -k =2,
otherwise

ork=n
ork=n-1
ork =n-2

h(k)x(n-k)=

k=-ao

From Fig. 2.13, note that

z(n

-k)-{~

Hence, the only nonzero terms in the sum are

y(n) = a"u(n) + 2a"-lu(n -1) + 3a"- 2u(n - 2)

The second way of evaluating the convolution is by using Eq. (2.26). The
evaluation, of course, produces the same result, but is simpler because of the

34

INTRODUCTION TO DIGITAL SIGNAL PROCESSING

{x(n)}

{h(n)}

..'I II. .

-3 -2 -1 0 1 2 3

.. 'I II...

-3 -2 -1

n=-1

n=-2

{x(k)}

{~~~~

(x(k)}

1II. ...

I 11. ... :

-4 -3 -2 -1

y(-2)=1
n=O

0 1 2 3

0 1 2 3 4 k

y(-I)=2

n=l

III
...

y(1)=2

y(0)=3

n=2

III .
III ..

(x(4)}

(h(2-k)}

-4 -3 -2 -1 0

(y(n)}.

I'l'I r

-4 -3 -2 -1 0

2 3 4 k

l'

y(2) ... 1

ftGURE loU
1bc CODvolution of two fiDite-duration sequences.

I' .

a.

1 2 3 4

DISCRETETIME DESCRIYnON OF SIGNALS AND SYSTEMS

3S

thIn)}
{x(n))

. rl I..
012

.... I

IIII I I

o1

II

{h(k)l
{x(n-k)}

.IJ

. . IIIIIII [

t"\N-C-

I I I

a:

s:

+I:

I ...

for n = 0

{x(n-k)l

. .1. .1. . . 1. ....._._.

'_00_ _. _

-2-1012

-\ 0 \ 2

for n = 1

for n = 2 ....

II I .

-10 1 2

{y(n)l

..... 1111 J1111 ,.;.

012

FIGURE 2-13
The convolution of a finite-duration sequence with an infinite-duration sequence.

finite duration of {x(n}}:

~

y(n) = .t ~
x(k) h(n - k) = l-O
~ (k + l)a"-tu(n - k)
__ oo

=a"u(n) + 2a"-lu(n -1) + 3a"-2u (n -

2)

Let us define a finite-duration sequence as one that contains only a finite

number of nonzero elements. When both {h(n)} and {x(n)} are finite-duration
sequences, then the convolution produces a finite-duration sequence, as
illustrated by Example 2.6. It is left as an exercise to show that if the input
sequence {x(n)} contains Nx elements and the unit-sample response sequence
{h(n)} contains NH elements, then the output sequence {y(n)} contains
Ny =N x + NH - 1 elements.

36

INfRODUcnON TO DIGITAL SIGNAL PROCESSING

A geometric series is a series in which consecutive elements differ by a

constant ratio. Such a series can be written in the form
x(n)

= r",

(2.31)

where r is a constant and Nl and N2 are any two numbers. In evaluating the
sum of products in the convolution, we commonly encounter sequences that
form such a series. The sum can then be easily computed by employing the
following two geometric sum formulas. To derive these formulas, let us
consider the following Taylor series expansion: Given a complex number c,
where Icl < 1, then

--=1+c+c2+c3 + ...
l-c

= ,,~
L c"

(2.32)

Applying this result in the reverse direction yields the infinite geometric sum
formula: Given a complex number c, where Icl < 1, then

L cit =1z

(2.33)
1- c
We can also compute the sum of a finite number of elements in a
geometric series. Let us consider the following sum
It-O

1+c+~+ +~-l=

N-l

,,0

C"

(2.34)

The sum of this finite-duration sequence can be expressed as the difference

between the sums of two infinite-duration sequences as

N-l
...
...
LC"=LC"-LC"
11==0

n-O

...

n-N

...

=Lc"-cNLc"
It-O

(2.35)

,,=0

Applying the infinite geometric sum formula to the two summations, we obtain
the desired finite geometric sum formula: Given a complex number c, then

N-l
l-cN
LC"=-It.O
1-c

(2.36)

While the infinite geometric sum formula requires that the magnitude of c
be strictly less than unity, the finite geometric sum is valid for any value of c.
EumpIe 2.8. Convolution of two lnfiIIite-dantion seqnences. Let
h(n) =a"u(n)
x(n) = b"u(n)

for all n,
for all n

DISCRETE-TIME DESCRIPTION OF SIGNALS AND SYSTEMS

h(n)=a"u(n) for all n

. .] rl~

.I(n)=/J'tu(n) for all n

11

to . . . .

.... IIIII

)0

-I 0 1 2

bb

-I 0 1 2 3

h(k)=t!u(k) for all k

. . ..'I rri

(.I(n-k)}
n=-1

It

iII: ..
JlIII1 r..

It

... 1JIII

-I 0

n=O

y(-I)=O
)0

It

y(O) = 1
)0

IlIlllH.

n=l

y(l)=a+b
)0

oI

n=2

)0

It

... I IIIII~H
1__ .

y(2)=a2+ab+b2
)0

012

It

y(n)

.I111I1I

)0

-10123

FIGURE Z-14
The convolution of two infinite-duration sequences.

as shown in Fig. 2.14. Then the output sequence values are equal to

y(n) =

L-

" akb"-k
aku(k)b"-ku(n - k) = L

k---

= b"

k-O

"

L (alb)"

k-O

The finite geometric sum formula can be applied to yield

r'

. ( )=b"I-(alb
n
t-(alb)

37

38

INTRODUcnON TO DIGITAL SIGNAL PROCESSING

2.6 STABILITY CRITERION FOR

DISCRETETIME SYSTEMS
For practical applications, we ace interested in designing only stable digit8I
filters, or those whose outputs do not become infinite. In this section, we
define a stability criterion and provide a necessary and sufficient condition for a
filter to meet this condition. In later chapters, we will describe alternate means
of determining whether or not a digital filter is stable.
Let us define a bounded sequence to be one in which the absolute value
of every element is less than some finite number M, i.e. {x(n}} is bounded if
Ix(n)1 <M for all n, where M is some finite number. Then, a system is said to
be stable when every possible bounded input sequence produces a bounded
output sequence. Since the output values of a digital filter should not exceed
the number range of the computer, it is appropriate to consider this bounded
output criterion.
STABll.1TY CRITERION. An LTI system is stable if, and only if, the stability
factor, denoted by S, and defined by

...

S=

L Ih(k)1
k=-CIO

(2.37)

is finite, i.e. S < 00.

The stability of a digital filter is expressed in terms of the absolute values
of its unit-sample response. We now show that this criterion is both sufficient
and necessary. To prove that it is sufficient, we show that if S is finite, then the
system is stable. Let {x(n}} be a bounded input sequence. We must show that
the output is bounded when S is finite. Taking the absolute value of the
output, expressed by the convolution equation, we have
ly(n)1 =

L~", h(k) x(n -

k)1

(2.38)

Recall that the magnitude of the sum of terms is less than or equal to the sum
of the magnitudes. Hence,
ly(n)1 ~

'"

Ih(k)x(n - k)1

(2.39)

k=-oo

Finally, since all the input values are bounded, say by M, we have for all n:

ly(n)1 ~ M

'"

Ih(k)1

k=-oo

or
ly(n)I~MS

(2.40)

Hence, since both M and S are finite, the output is also bounded.
To show that it is necessary for S to be finite for the filter to be stable, we
present a counter-example. We need to finq a bounded input sequence that

DISCRETETIME DESCRIPfION OF SIGNALS AND SYSTEMS

39

produces an unbounded output. For our input sequence, let us choose

x(n)

= {~( -n )/lh( -n)1

for all n for which h( -n)

otherwise

'* 0

(2.41)

This sequence consists of plus and minus ones and is obviously bounded. Let
us now consider the output value for n = 0, or y(O). From the convolution
equation, we have

...

y(O) =

2:

...

h(k)x( -k) =

k--~

h2(k)/lh(k)1

k=-~

...

2:

(2.42)

Ih(k)I=S

k==-oc

Hence, for y(O) to be bounded, it is necessary for S to be finite. We have thus

proved the necessary and sufficient conditions. The following examples
illustrate the application of the stability criterion.
Eumple %.9. StabiHty of the Ceaeralized J-sample avenger. Let us generalize
the 3-sample averager by letting the coefficients be arbitrary constants b -1> b o and
b l . The output is then equal to
y(n)

= b_Ix(n + 1) + boX(n) + blx(n -1)

To test whether this filter is stable, we compute the value of S:

. --

s= ~

Ih(n)1 = Ib_11 + Ibol + Ibll

Hence, S is finite for finite values of b_ .. b o and b l .

Example 2.9 can be easily extended to show that any filter with a
unit-sample response having a finite number of nonzero elements is always
stable.
Eumple 1.10. Stability of first-onler reamive &her. Let
y(n)

=ay(n -1) +x(n)

Recall that h(n) =a"u(n), for all n. Checking for stability, we must calculate

s = ,,_-OD
~ Ih(n)1 = ~ lal"
11-0

It is obvious that S is unbounded for lal2: 1, since then each term in the series is
2:1. For lal < 1, we can aIfply the infinite geometric sum formula, to find

1
S=--

1-a

forlal<1

Since S is finite for lal < 1, the system is stable for lal < 1.

40

INTRODUCTION TO DIGITAL SIGNAL PROCESSING

Example Z.ll. Stability of seoond-order filter. Let

yen) = 2r cos(wo) yen - 1) -,z yen - 2) + x(n) - r cos(wo) x(n -1)

To test the stability of this filter, we must first determine its unit-sample response
and then compute the value of S. To determine {hen)}, we assume zero initial
conditions, or yen) = 0, for n <0 and let the input sequence be {den)}. We then
get the following output values.

yeO) = 1
y(l)

=2r cos(wo) -

y(2) =

r cos(wo)

=r cos(wo)

2r'- cos ( wo) - ,z = ,z(2 co~( wo) - 1) = r2 cos(2wo)

2

y(n) = r" cos(nwo)

Hence,

h(n) = ", cos(nwo) u(n)

Checking for stability, we have
-

00

~ Ih(II)1 = ~ I'" cos(lIwo)1 s ~o Irl" = l-Irl

for Irl < 1

Hence, the filter is stable for Irl < 1. The inequality of the sums in the
above stability calCUlation comes about because the magnitude of the cosine
function is less than or equal to unity.

2.7 CAUSALn:Y CRITERION FOR

DISCRETE-11ME SYSTEMS
A system is said to be causal when the response of the system does not precede
the excitation; All physical systems are causal in that they do not react until a
stimulus is applied. An LTI system is causal if, and only if,
hen) =0

for n <0

(2.43)

:his is true because, by definition, {hen)} is the response of a system to a

unit-sample sequence, whose nonzero element occurs at n = O. For a causal
system, the output cannot have nonzero values for n < 0, the time before the
nonzero-valued element is applied.
Causal filters are typically employed in applications in which the data
samples are processed as they are received. Qearly, in this case we cannot
have access to future sample values. In this book, we will not restrict ourselves
to causal digital filters. When filters are used to process data for which the
entire data record is available at the start of processing, there is no restriction
against looking ahead to future values. One example of such data is that
obtained from images. As described in the previous chapter, this look-ahead
feature is implemented by the advance, or negative delay, component of the
digital filter. For noncausal filters, {hen)} has nonzero elements for n <0, as
illustrated with the 3-sample averager in Example 2.3. Since the unit-sample

DISCRETE-TIME DESCRIPTION OF SIGNALS AND SYSTEMS

41

response of the first-order recursive filter in Example 2-4 obeys the condition
given in Eq. (2.43), it is a causal filter.

2.8 LINEAR CONSTANT-COEFFICIENT

DIFFERENCE EQUATIONS
In this section, we formalize the notation that we have been using above for
describing the time-domain behavior of digital filters. In general, the output of
a finite-order LTI system at time n can be expressed as a linear combination of
the inputs and outputs in the following form:
M

Np

y(n) = ~ a/cy(n-k)+ ~ b/cx(n-k)

k=l

(2.44)

k=-NF

In words, the current output y(n) is equal to the sum of past outputs, from
y(n -1) to y(n - M), which are scaled by the delay-dependent feedback
coefficients ak, plus the sum of future, present and past inputs, which are
scaled by delay-dependent feedforward coefficients bk The output value at
time n is determined by the input values from N F samples in the future
x(n +NF), through the current value x(n), to the sample N p time units in the
past x(n - Np ). For a causal filter NF sO. When NF> 0, then the future values
of the input determine the current output value, and the filter is noncausal.
EumpIe 2.12. Difference equation of a simple filter. We can combine the
equations of the nonrecursive 3-sample averager and the first-order recursive
filter, to get
y(n) = ay(n -1) + lx(n + 1) + lx(n) + lx(n -1)
This filter has both a recursive and a nonrecursive part and is also noncausal.

The difference equation is useful not only for analyzing digital filters but
can be employed for implementing them. A digital filter can be implemented
from the general difference equation by performing the following four steps, as
shown in Fig. 2.15.

Step 1. Two points are drawn, one corresponding to the current input
x(n) and the other corresponding to the current output y(n).
Step 2. To gain access to the past values of the input and output
sequ~nces, delay elements are connected to the output and input points.
Future values of the input sequence are obtained by connecting advances to
the input point.
Step 3. Multipliers are connected to the outputs of the delay elements to
produce the required products. Multipliers with gains of unity need not be
included, but replaced by a simple wire.
Step 4. Finally, the digital filter is accomplished by connecting the
outputs of the multipliers to two adders. The first generates the sum of

42

INTRODUcnON TO DIGITAL SIGNAL PROCESSING

STEP I
x(n)_

y(n)'

STEP 2

..+N,)--+

x(n+2)~

,..---.. y(n)
x(n)---~

I
I
I

~
--..T~y(n-M}
STEP 3

y(n)
x(n)--+----i

ftGURE 2-15
Steps to produce a digital filter structure from the general difference equation. This filter structure

will be implemented in the computer programs.

DISCRETE-TIME DESCRIPTION OF SIGNALS AND SYSTEMS

43

STEP 4

X(II)

_-.------1 >-..--,~

}-----_--oY(II)

products of the inputs and the feedforward coefficients, and the other that of
the outputs and the feedback coefficients.
The following example illustrates the application of these four steps.

_tal

Example 1.13. lmpIemeatatioD of the seeoad-order

filter fro. the
dlrereace eqaatioa. The filter implementation procedure from the difference
equation given in Example 2.11 is shown in Fig. 2.16.

X(II)_-.......----~

}------.---y(II)

FIGURE 2-16
Digital filter structure that implements the difference equation given by
y(n) = 2r cos(ClJo) yen -1)

_,2 yen -

2) + x (n) - , cos(wo)x(n -1)

44

INTRODUCTION TO DIGITAL SIGNAL PROCESSING

x(n) o---+----i

.>--30(

J..------__-- y(n)

ftGURE 1-17
General digital filter structure implemented in the computer programs.

The general digital filter structure shown in Fig. 2.17 is convenient for
describing the logic behind the computer programs described below. However,
it is not the most efficient structure in terms of the number of delays that are
used. In a later chapter, we will consider other filter structures that have a
smaller number of delays but for which the corresponding computer programs
are more complicated. Conceptually, however, these more memory-efficient
programs perform the same types of operations as the simple programs
described below.

2.9 SUGGESTED PROGRAMMING

STYLE
Befqre describing the programs, we provide the motivation for the particular
programming style suggested below. The goal of the programs is to yield the
greatest intuitive understanding of the digital filter operation with the least
programming effort. Techniques to improve the efficiency of these programs
are suggested throughout the text. The logical steps in the operation of a
digital filter are best illustrated by writing simple subprograms that perform the
primitive tasks that are most commonly employed in signal processing. These
simple subprograms are then combined to form the digital filter program.
To provide a FORTRAN example of the recommended programming
style, one common task is to initialize an array to zero. Even though this can
be done quite easily in a simple loop in the main proram, the programming

DlSCRETETIME DESCRIPTION OF SIGNALS AND SYSTEMS

4S

style suggested above would use the following subroutine:

SUBROUTINE NULL(X,NX)
DIMENSION X(NX)
DO 1 I=1,NX
X(I)=O.
RETURN
END

An array B(N), for 1:s N:s NB, can be zeroed in the main program with the
following statement:

CALL NULL(B,NB)
Two additional subroutines, described below, compute the sum of
products of two arrays and implement a shift register. Use of such simple
subroutines will simplify the main program by making it more readable, which
is important for debugging and trouble shooting. This programming style also
allows the student to reason in terms of these primitive tasks.
In trying to explain the programs in a manner that is not specific to a
particular programming language, a problem immediately arises with aI'tay
indexing. The various computer languages used for signal processing can be
divided into two classes: those that allow arrays to be indexed starting with
zero, such as BASIC and C, and those that do not, such as FORTRAN and
PASCAL. In the first set, x(O) is the first element of an array, in the second.
the first element is x(l). Since all languages can start indexing from 1, we will
use this convention in the programs below. For computer languages that allow
the array index to start from zero, a slight modification in the indexing of the
programs should be made. The advantage of starting the array with x(O) is that
it then conforms to the notation used for sequen~es in the analytic sections of
this book, i.e. x(n), for 0 :s n :s N - 1.

2.10 WRITING A DIGITAL FILTER

PROGRAM
In this section, we describe the cOJ1lputer programs to implement a digital filter
from a difference equation. These programs will be used in the following
chapters for analyzing and implementing digital filters to meet some desired
specification. To implement a digital filter as a flexible program, it is most
convenient to write separate subprograms for the nonrecursive part, to be
called NONREC, and the recursive part, to be called REC. The general filter
can then be implemented by first using NONREC followed by REC, as
described below.
First, it is useful to establish the notation for the arrays to be used in the
programs. We are typically given an input array of real values, denoted by

46

INTRODUcnON TO DIGITAL SIGNAL PROCESSING

X(N), for 1 s; N s; NX, where NX is a known value. This array represents the
input data to be processed. In this book, the most common input will be the
unit-sample sequence. We are usually asked to design a digital filter that is
specified by an array of feedback coefficients ACOEF(N), for 1 s; N s; NA, to
be used in REC and an array of feedforward coefficients BCOEF(N), for
1 s; N s; NB, to be used in NONREC.
Since we are implementing nonrecursive and recursive filters separately,
it is helpful to define two output arrays. The output of a nonrecursive filter is
denoted by YN (N) and will have a duration NY = NX + NB - 1. This duration
follows from the convolution of two finite-duration sequences. The output
array of the recursive filter is denoted by YR(N) and has a duration NY. For
the recursive filter, however, the value of NY must be specified. Because of
the feedback loop, the output of an ideal recursive filter contains an infinite
number of non-zero elements. However, for a stable filter, the output values
eventually approach zero when the input has a finite duration. Hence, there
will be only a finite number of output elements whose amplitudes are
significantly greater than zero. For all the filters in this book, NY will be less
than or equal to 128.
Since most programming languages require that the size of all arrays be
specified at the beginning of the program, a length of 128 floating-point, or
Teal, elements for all arrays will be adequate for most of the projects in this
book. In the descriptions of the variables below, Teal is used to specify a
floating-point variable.
We first consider the implementation of a causal filter, for which NF=O
in Eq. (2.44), and then describe the simple modification that can be added to
implement a noncausal filter.
PROGRAMMING A CAUSAL FILTER. The general causal filter is implemented by combining two causal components: a strictly nonTecuTsive causal
part with a strictly TecuTsive part.
The general difference equation (2.44) of a causal nonTecuTsive filter, for
which NF = 0 and ak = 0 for all k reduces to
Np

yen) =

bkx(n - k)

(2.45)

k=O

The feedforward coefficients bk for 0 s; k s; Nr are stored in BCOEF(N) for

1 s; N s; NB (NB = Np + 1). The value of bo is stored in BCOEF(l), b l in
BCOEF(2), etc. Notice that Eq. (2.45) has the form of the following
convolution
Np

y(n) =

L h(k)x(n -k)

(2.46)

k=O

By computing the unit-sample response of the nonrecursive filter, we find that

it has a finite duration and is equal to the feedforward coefficient sequence, or
h(k) = bk

for - NF s; k

s; Nr

(2.47)

DISCRETETIME DESCRIPTION OF SIGNALS AND SYSTEMS

47

Hence, a nonrecursive filter can be implemented directly from its unit-sample

response. For the causal filter considered here, Np = O.
The difference equation for a strictly recursive filter, whic~ employs the
past outputs and only the current input value, for which bk = 0 for k =1= 0,
reduces to
y(n) =

L aky(n-k)+x(n)

(2.48)

k-l

The feedback coefficients ak for 1 s k s M will be stored in ACOEF(N) for

1 s N s NA, where NA = M. The value of a 1 is stored in ACOEF(1), az in
ACOEF(2), etc.
In both Eqs. (2.4S) and (2.48), a sum-of-products calculation employing
an array of coefficients and an array of data is required. We implement this
primitive task as a separate subprogram, called SUMP (sum of products), that
can be used in both NONREC and REC. To simplify the programming logic,
we will use two auxiliary buffer arrays. The first, denoted XBUF(N), for
1 s N s NB, stores the current and past NB values of {X(N)}. Only NB values
are needed since that is the number of coefficients in BCOEF(N). The second
buffer array, denoted YBUF(N), for 1 s N:5 NA, stores the past NA values of
{Y(N)}. A second subprogram, called SHFTRG (shift register), acts as a shift
register to maintain the proper time sequence of values in these arrays.

SUMP Subprogram. The SUMP subroutine computes the sum of products of

two arrays, a coefficient array, denoted COEF(N), for 1sN:5NC, and a
buffer array BUF(N), of the same size, to produce the output value SUM. The
call to the subprogram is then

SUMP(COEF,BUF,NC,SUM)
where
COEF
BUF
NC
SUM

(real input array) is an array of coefficients;

(real input array) is an array of data values;
(integer input) is the duration of the COEF array;
(real output) is the sum of products equal to
NC

SUM = L COEF(I)

* BUF(I)

(2.49)

1-1

as shown in Fig. 2.18. When SUMP is used in NONREC, BUF is equal to

XBUF, COEF corresponds to BCOEF and NC is equal to NB (= Np + 1).
When used in REC, BUF is equal to YBUF, COEF corresponds to ACOEF
and NC is equal to NA ( = M).

SIIFI'RG Subroutine. To update the buffer arrays, we implement a shift

register routine, called SHFTRG, that shifts the contents of the buffer array by
one element at each iteration of the filter operation. The shift is performed

48

INTRODUCi10N TO DIGITAL SIGNAL PROCESSING

COEFFICIENT ARRAY

BUFFER ARRAY

COEF(I)

BUF(I)

COEF(2)

BUF(2)

COEF (3)

BUF(3)

COEF(NC)

BUF(NC)

nGURE 218
Operation performed by the SUMP
routine.

such that the oldest element in the buffer is lost and the current value of the
time sequence is put into the first element of the buffer. That is, at time L,
X(L) is put into XBUF(1) and only the current and past NB -1 values are
retained in the buffer.
In general, given a value VAL and buffer array BUF(N), for 1::5 N ::5
NBUF, we want to write the subroutine called
SHFTRG(VAL,BUF,NBUF)
where
VAL
(real input) is the value to be inserted into the buffer;
BUF
(real input and output array) is the buffer array;
NBUF
(integer input) is the duration of the buffer array.
This subroutine shifts BUF(K) into BUF(K + 1), for 1::5 K::5 NBUF, and puts
VAL into BUF(l), as shown in Fig. 2.19. When used in NONREC, the buffer
BUF is equal to XBUF, NBUF is equal to NB, and at time L, VAL is equal to
X(L). When used in REC, BUF is equal to YBUF, NBUF is equal to NA and
at time L, VAL is equal to YR(L). We now describe how to use these
subroutines to implement the nonrecursive filter subprogram NONREC and
the recursive filter filter subprogram REC.
NONREC Subprogram. For the nonrecursive filter, the input array X(N), for
1::5 N::5 NX and the feedforward coefficient array BCOEF(N), for 1::5 N::5 NB
BEFORE SHFfRG
CALL

AFfER SHFfRG
CALL

VAL
BUF(I)

BUF(I)

BUF(2)

BUF(2)

BUF(NC-I)

BUF(NC-I)

BUF(NC)

BUF(NC)

nGURE 2-19
Operation performed by the SHFrRG
routine.

DISCRETE-TIME DESCRIPTION OF SIGNALS AND SYSTEMS

49

(NB = Np + 1), are specified and the output array YN(N), for 1:5 N:5 NX +
NB - 1, is to be generated_ The nonrecursive filter subroutine is denoted by

NONREC(X,NX,BCOEF,NB,YN)
where
X
NX
BCOEF
NB
YN

(real input array) is the input array to the nonrecursive filter;

(integer input) is the duration of the input array;
(real input array) is the feedforward coefficient array;
(integer input) is the duration of the coefficient array;
(real output array) is the output array of duration NY = NX +
NB-L
The procedure to be performed within NONREC is indicated in the logic
diagram shown in Fig_ 2_20. The NONREC routine is accomplished using the
three basic steps described below.
Step 1. The buffer array XBUF(N), for 1:5 N:5 NB, is defined and
initialized to zero.
Step 2. In a loop with index I, for I equal from one to NX, VAL is set
equal to the Ith element of the input, or X(I). VAL is then entered into the
shift register XBUF. The sum of products of XBUF and BCOEF are
computed to produce SUM. The Ith value of the output sequence for the
nonrecursive filter, YN(I), is then set equal to SUM.
Step 3. When I> NX, all the elements of X(N) have been used. VAL is
then set to zero and then entered into the buffer XBUF. The sum of products
is computed and the result SUM becomes the value of the output YN(I). This
loop is executed until no more elements of the input sequence X(N) remain in
the buffer, or until I > NX + NB - 1.

After this routine is completed, the nonrecursive filter operation has been
accomplished and the output values will reside in the YN array. The
combination of the SHFTRG and SUMP routines results in a convolution
operation.
REC Subroutine. For the recursive filter, the input array X(N), for 1:5 N:5

NX, and the feedback coefficient array ACOEF(N), for 1:5 N:5 NA (NA =
M), are specified and the output array YR(N) , for 1:5 N:5 NY, is to be
generated. Unlike the nonrecursive filter, the duration of a recursive filter can
be infinite because of the feedback. Hence, the value of NY must be
sufficiently large so that the important features of the output are observed. In
the case of determining the unit-sample response, the value of NY must be
chosen large enough such that the values of h(n) for n > NY are approximately
zero. In producing a plot of the sequence, it is sufficient that the omitted
magnitudes be smaller than the resolution of the plotting device. Typical
values for NY used "in the projects will be 55 and will not exceed 128.

50

INTRODUcnON TO DIGITAL SIGNAL PROCESSING

FIGURE 2-20
Logic diagram for the NONREC filter routine. The values for the input array X(N), the size of the
array NX, the feed-forward coefficient array BCOEF(N), and the size of the coefficient array NB
are given; the output array of the nonrecursive filter YN(N) is to be computed.

DISCRETE-TIME DESCRImON OF SIGNALS AND SYSTEMS

51

The recursive filter is implemented by the subroutine

REC(X,NX,ACOEF,NA,YR,NY)
where
X
(real input array) is the input array to the recursive filter;
NX
(integer input) is the duration of the input array;
ACOEF
(real input array) is the feedback coefficient array;
NA
(integer input) is the duration of the coefficient array;
YR
(real output array) is the output array;
NY
(integer input) is the duration of the output array.
The logic diagram for the REC routine is shown in Fig. 2.21. Within REC, the
following steps are performed:
Step 1. The buffer array YBUF(N), for 1:5 N:5 NA, is defined and
initialized to zero.
Step 2. In a loop with index I, for 1:5 1:5 NX, the output at time I,
YR(I), is equal to the sum of the input X(I) and SUM, the product of the past
outputs stored in YBUF and ACOEF. For 1= 1, SUM is equal to zero. VAL
is then set equal to the output YR(I) and is sequenced into YBUF.
Step 3. When I > NX, all the elements of X(N) have been used. At
that point the output YR(I) is then determined from only the past output
values. The recursive filter operation is accomplished after the desired number
of output points NY have been computed.
GENERAL FlLTER IMPLEMENfATION. The general filter can be imple-

mented by calling the NONREC followed by the REC routines, as shown in

Fig. 2.22. The input sequence X(N), for 1:5 N:5 NX, to NONREC produces
the output sequence YN(N) , for 1:5 N :5 NX + NB - 1. This sequence then
becomes the input sequence to the REC routine. The REC routin~ will then
produce the output sequence YR(N) , for 1:5 N :5 NY. Usually NY will be
greater than NX + NB - 1. The array YR is then the output of the general
filter.
If the REC routine is called first and then followed by NONREC, care
must be taken in choosing the value of NY, the duration of the output
sequence YR produced by REC. This sequence then becomes the input to
NONREC. The output of NONREC, YN, then contains NY + NB - 1
elements, a duration that must not exceed the dimension of YN.
NONCAUSAL FILTER IMPLEMENfAnON. By definition, a filter is noncausal

when future values of the input sequence are used in the computation of tb,e
current output value. Noncausal filters can be identified easily from the
difference equation by Np being greater than zero. Since arrays with negative
indices are not generally allowed in programming languages, the simplest way
to implement a noacausal filter is to change the difference equation in such a

52

INTRODUcnON TO DIGITAL SIGNAL PROCESSING

~{

nGURE 2-21
Logic diagram for the REC filter
routine. The values for the input
array X(N), the size of the array NX,
the feedback coefficient array
ACOEF(N), and the size of the
coefficient array NA are given; the
output array of the recursive filter
YR(N) is to be computed.

way as to produce a causal filter and then shift the final output array into the
future by N p index units, as described below. Using this procedure, the causal

filter can then be implemented in the manner described above, and only an
additional shift routine is needed, to obtain the noncausal result.
The noncausal property of the filter is due to the NF feedforward
coefficients having negative index values in the nonrecursive part of the filter.
Hence, a noncausal filter can be made causal by shifting the coefficient

DISCRETETIME DESCRIPTION OF SIGNALS AND SYSTEMS

YN(N)
lo;;No;;NX+NB)

NONREC
X(N)
for lo;;No;;N

:x"

53

B.E.C
YR(N)

ACOEF(N) - 1o;;No;;NY
lo;;No;;M

BCOEF(N)
lo;;No;;NB

FIGURE 2-U
Combining NONREC and REC to form the general filter having both nonrecursive and recursive
parts.

sequence {b k } to produce {ba, where

as shown in Fig. 2.23. This new set of coefficients then generates the causal
output, denoted by {y'(n)}, whose elements are equal to
M

y'(n)=

Np+Np

L aky'(n-k)+ L
k-l

b;'x(n-k)

(2.50)

k~O

NON CAUSAL FILTER

CAUSAL FILTER

Noncausal
y'(n)

part

zen) c---~---l

~~-~

(y(n)}

,... rII"r

(ytn)} =(y(n.NF)}

, .. 1 I lour+ ~n

o
Ji1GURE 2-23
Transforming a noncausal filter into a causal filter.

NF

Np

54

INTRODUCfION TO DIGITAL SIGNAL PROCESSING

CAUSAL FILTER OUTPUT

{/(n))

oorIIIlhr. ..
0

~7
I

NON CAUSAL (y(n)}

FILTER OUTPUT

1 2 NF NF+Np

Jr1 rT

_ .......~T';-L-L-L-!-'L.JWW'L..o...........o--_ _ _ _ _.".:>

- NF

Np

FIGURE 2-24
The noncausal filter output is
obtained by shifting the output
of the causal filter to the left.

This causal tiltf!r can be implemented by using the routines described above.
The values of {y'(n)} are exactly the same as those of {y(n)}, but only shifted
by NF index values. To obtain the noncausal output sequence {y(n)}, a time
shift to the future is required, i.e.
y(n) = y'(n + NF ),

for -NF ::; n::; NY - N F

as shown in Fig. 2.24. This time shift is performed by the subprogram

SHIFT(Y,NY,NF)
where
y

(real input and output array) is the array to be shifted;

(integer input) in the duration of the array;
(integer input) is the number of elements that the input array is to
NF
be shifted to produce the output.
After the shift is performed, the first NF samples of the input array are lost
and the last NF samples of the output array are set to zero, as shown in Fig.
2.25.
NY

AFTERSHlFT
CALL

BEFORE SHIFT
CALL
Y(l)

Y(l)

Y(2)

Y(2)

Y(NF )

Y(Ny-Nr-l)

Y(NF + 1)
Y(NF +2)

Y(Ny-NF)
Y(Ny-NF+l)

Y(Ny-l)

Y(Ny -1)

Y(Ny)

Y(Ny)
0

FIGURE l-2S
Operation performed by the SHIFf
routine.

DISCRETE-TIME DESCRIPTION OF SIGNALS AND SYSTEMS

55

In programming this shift, we still have a problem, because the shifted

array may require negative index values_ To avoid this problem and still
observe the noncausal components, if any, in the filter unit-sample response,
the input to the filter can be the delayed unit-sample sequence, whose nonzero
component occurs at N = NF Using this as the input sequence, the noncausal
components of the filter can then be observed in the output sequence YN(l)
through YN(NF -1).

2.11 SUMMARY
This chapter has presented the time-domain analysis techniques for describing
linear time-invariant digital filters and discrete-time systems, from basic
definitions to programming techniques. A particularly useful sequence that
describes the behavior of a digital filter is its unit-sample response, since from
it we can determine whether a filter is stable and whether it is causal. Further,
for any input sequence to a digital filter, we can determine the output sequence
by performing convolution with the unit-sample response.
An alternate time-domain description of a digital filter is the linear
difference equation. Yhe unit-sample response can be determined from the
difference equation by simply making the input equal to the unit-sample
sequence and computing the output sequence. Of practical importance, it was
found to be a rather straightforward task to implement a digital filter structure
from the difference equation. If the unit-sample response has a finite duration,
a nonrecursive filter can be implemented by identifying the feedforward
coefficients with the corresponding elements of the unit-sample response.
The steps for implementing a digital filter as a computer program were
presented. By writing simple routines to perform a sum-of-products computation and a shift-register operation, the implementation of nonrecursive and
recursive filters becomes conceptually easy to grasp. The general filter
structure containing both nonrecursive and recursive parts is implemente.d by
employing both routines. Noncausal filters can be implemented by employing
the causal filter routines and applying a subprogram that shifts the output
sequence "into the future".

PROBLEMS
2.1. Express the following in complex vector notation:
(a) 1 + e- j8
(b) 1- a e- j8
(c) 1 + 2e- j8 + e-j28
(d) 1 + ja
1 + jb

1 + a eJ8
(e) 1 + b eJ8
(f) (a

+ jb )112

56

INfRODUCfION TO DIGITAL SIGNAL PROCESSING

2.2. Express the following in real plus imaginary notation:

(a) 2cos(0)e- Jw
(b) _I_
1+ ja
1
(c) 1 + a ej8

(d) 1 + e

J; + e J28

(e) (20 + 16 cos 0)112 exp[-j arctan( -sin 0 )]

2+cos~

2.3. Use the Euler identities to derive the formula for the cosine of the sum of two
angles:
R) cos( a + P) + cos( 11'- P)
cos( 11'+ P =
2
2.4. Demonstrate that the first-order recursive filter is linear. (Hint: Express G{x(n)}
as an infinite sum of scaled terms containing current and past values of x(n).
2.5. Evaluate the following summations and, give ,the result in complex vector
notation:
for Irl < 1
n=O

N-)

(b)

L e-

jwn

n=O
N-)

(c)

L cos(wn)

n=O
N-I

(d)

L an e-

jwn

2.6. This problem considers the case of determining the unit-sample response of a
filter having nonzero initial conditions. This condition can occur when the delay
element contents are not initialized to zero.
(a) The 3-sample averager. Determine {h(n,O)} when the content of the delay
element is equal to 2 and that of the advance is equal to 1 at n = O.
(b) The first-order recursive filter. Determine {h(n,O)} when the content of the
delay element is equal to 3 at n = O.
2.7. Given two finite-duration causal sequences, {x(n)} containing Nx elements and
{h(n)} containing NH elements, show that {y(n)} = {x(n)} * {h(n)} produces a
causal, finite-duration sequence containing Ny = Nx + NH -1 elements.
2.S. Given the following two sequences, determine an analytic expression for
{y(n)} = {x(n)} * {h(n)}. Let
x(n) = {~

forOsnsN-1
otherwise

{~

for OSnsN-1
otherwise

and
h(n) =

DISCRETE-TIME DESCRIPTION OF SIGNALS AND SYSTEMS

2.9.
2.10.

2.11.

2.U.

2.13.

2.14.

2.15.

2.16.
2.17.
2.18.

2.19.

57

Hint: First try N = 2 and then N = 3, to obtain the general analytic form of the
sequence.
Prove that any realizable filter having a finite number of elements in its
unit-sample response is always stable.
For the following difference equations, determine and sketch the unit-sample
response sequence and implement the digital filter structure:
(a) y(n)=x(n)-x(n-N)
(b) y(n)=ay(n-1)+x(n)+x(n-1)
Let a = 0.9.
For the digital filter structures shown in Fig. n.ll, determine the difference
equation, compute and sketch the unit-sample response.
The steady-state condition of a discrete-time system is defined as that in which the
output value does not change with time, i.e. y(n + 1) = y(n) for n very large.
Determine this steady-state value when the input sequence is a constant-valued
step sequence, x(n) = Ku(n), for the following system:

y(n) = ay(n -1) +x(n) + x(n -1)

Let a = 0.9.
Implement the digital filter that generates cos( won), for n ~ O. (Hint: see
Example 2.11).
Demonstrate that the sequence {x(n)} = {sin(wn)} produces identical values
when w = Wo and w = Wo + 2mJr, where m is an integer. Sketch the results for a
particular value of wo, with m = 0 and 1.
Demonstrate that the convolution operation is commutative by performing the
convolution given in Example 2.8 in the reverse order to obtain the same result.
Restructure the logic diagram of the nonrecursive filter such that only one loop is
required.
Restructure the logic diagram of the recursive filter such that only one loop is
required.
For each of the following transformations, determine whether the system is
stable, causal, linear and shift-invariant:
(a) Gt{x(n)} = ax(n - no) + bx(n - nt), for constants a, b, no and nl
(b) G2 {x(n)} =x~n)x(n - no)
(c) G3 {x(n)} =x(2n)
That a system is unstable does not imply that the output at time n, or y(n),
approaches infinity for every input sequence. Consider the system whose
unit-sample response is given by

h(n)={~

for n ~O
otherwise

(a) Show that this system is unstable.

(b) Let the input sequence be given by
x(n) = eJwn

for all n

Using the convolution equation, determine the values of w for which y(n) is
finite, and the values of w for which y(n) is infinite.

58

INTRODUcnON TO DIGITAL SIGNAL PROCESSING

(a)

x(n) - ___---~

~---___oy(n)

(b)

x(n)c

(c)

x(n)

y(n)

y(n)

(d)

y(n)

x(n)

(e)

x(n)

y(n)

OGURE Pl-ll
Digital filter structures.

2.20. Find the output sequence produced for the system having the unit-sample
response given by

h(n)=H

forn=O
forn=4
otherwise

DISCRETETlME DESCRIPTION OF SIGNALS AND SYSTEMS

S9

when the input sequence is given by

forn<:!O
otherwise
2.21. A common input to an LTI system used in control situations is the unit-step
sequence, {u(n)}. The response of the system is the unit-step response, denoted
by {g(n)}. Determine the unit-step response of the first-order recursive system,
whose unit-sample response is given by h(n) = a"u(n).
x(n) =

{~

COMPUTER PROJECTS
Preface
Most of the projects contain the following parts:
Object, that is a concise statement of the purpose of the project.
Analytic results, which describe the required analysis to be performed for
later comparison with the computer results. In some cases, references are
made to a problem at the end of the chapter.
Computer results, which describe the programs to be written and the
computations to be performed.

Definitions for the following terms will be helpful for performing the projects.

AnalyticaUy determine means to apply the appropriate mathematical

formulas and evaluate the expression. All of the summations in the analytic
expressions should be evaluated and complex quantities should be expressed in
complex vector notation. The analysis is usually not difficult when requested. If
the student finds that the procedure is becoming very complicated, it means
that either a mistake has been made, or that a simpler approach can be taken.
In this case, the programming should be performed first to provide the answer
towards which the analysis should lead.
Sketch means to provide a graphical interpretation of the analytic results.
If the results were not analytically determined, the graphical interpretation is
to indicate a qualitative understanding of the operation being discussed. In all
cases, the time- or frequency-domain scales and landmarks are to be indicated
on the sketches. In many cases, only relative amplitudes are necessary.
Compute means that the procedure is to be performed numerically by the
computer and compared with the analytic results.
Plot means that the computed data are to be presented in graphical form
and compared with the sketches obtained from analysis. The graphs produced
by the computer should be labeled by hand to indicate the correct time scale
(n) or frequency scale (co). The amplitudes and abscissa locations of features

60

INTRODUcnON TO DIGITAL SIGNAL PROCESSING

of interest in the plots, such as peaks, valleys or discontinuities, should be

noted and related to the sketches.
The project report should be submitted with the material for each project
grouped together. The analysis should be presented first, followed by the
corresponding computer plot for easy comparison. Lists of numerical printouts
should not be included. The programs used in the project should be appended
to the end of the report, to facilitate evaluation of the programming style and
for easy reference in the future.
2.1. Sinusoidal sequences
Object. Demonstrate that the values of sinusoidal sequence are unique for
frequencies only in the range -1r:S W < 1r.
Analytic results. Demonstrate that {x(n)} = {sin(wn)} produces identical values
when w = Wo and w = Wo + 2m1r, where m is any integer and - n :s Wo < n. Sketch
the results for a particular value of wo, with m = 0 and 1.

Computer results. Using the program TESTO.FOR in Appendix A, enter the two
values of frequency used above to produce the sinusoidal sequences. The program
asks for the frequency in Hertz per sample interval F, where F = w /2n. The range
of F that produces unique values is then -0.5:s F < 0.5.
Use the PLOT.FOR subroutine given in Appendix A to plot the data on the
terminal and on the printer. Since 55 lines fit on one piece of paper, limit the
number of points plotted to 55 whenever possible.
2.2. Exponential sequences
Object. Generate an exponential sequence and to determine the value number
range of the computer.

(1) Real-valued exponential

Analytic results. Determine and sketch the exponential sequence
x(n) = r"

for 0 :s n :s 8.

for some positive value of 0.9 < r < 1. Repeat for some negative value of
-l<r< -0.9.
Computer results
(a) Write a program to accept a real-valued constant r and compute 55 elements
of the exponential sequence. Plot the sequences for the two values of r used
above. Compare your results with analytically determined sequences and
sketches.
(b) Determine an approximate range for the largest floating-point number that
can be expressed on your computer by accepting increasingly larger values of r
until an overflow eITor is encountered in computing the sequence. Determine
the smallest magnitude number that is still nonzero by accepting decreasingly
smaller values of r, until either an underflow eITor occurs, or the sequence
values become zero.

DISCRETETIME DESCRIPTION OF SIGNALS AND SYSTEMS

61

(2) CompLex-vaLued exponential

Generate and plot 55 elements of the real and imaginary parts for the complex
exponential sequence given by
x(n) = (r e JWo)"

for Wo = :rIB and for the two values of r used above.

2.3. Determination of the effective duration of discrete-time sequence
Object. Determine the effective duration of an infinite duration unit-sample
response sequence.
For the infinite-duration exponential sequences, the magnitude of the
elements typically decreases as rn. Since the plotting routines have a resolution of
approximately 1 percent, we can define an effective duration of the sequence, as
the number of sequence values that are greater than 1 percent of the maximum
value of the sequence. Beyond this effective duration, all the elements can be
considered as being zero-valued.
Analytic results. For the real-valued exponential sequence x(n) = rnu(n),
determine an analytic expression for the effective duration of the sequence as a
function of r.
Computer results. Verify this result by plotting several representative sequences.
2.4. Verification of the digital filter programs
Object. Write and verify the subroutines to implement nonrecursive and recursive
filters.
Write the subroutines NONREC and REC as described in the text. To
implement the noncausal feature, write the subroutine SHIFT. To verify the filter
operation, let the input array be a delayed unit-sample sequence, whose nonzero
element occurs at N = 6, i.e., X(6) = 1, X(N) = 0, for 1 =s N =s NX, N 6. With
this definition, the sixth element of the output sequence will correspond to n =
and the noncausal components will occur in Y(l) through Y(5).

For each of the filters below, perform the following:

(i) draw the filter structure;
(ii) determine an analytic expression for the unit-sample response sequence and
sketch the sequence;
(iii) compute and plot the unit-sample response;
(iv) compare the analytic and computed results, indicating the agreement or
explaining any differences.
Computer results
(1) Causal nonrecursive fiLters. Implement the following difference equations as
causal digital filters. Let NX = 20.
(a) y(n) =x(n) + 2x(n -1) + 3x(n - 2)
(b) y(n)=x(n)-x(n-lO)
(2) NoncausaL nonrecursive fiLter. Implement the 3-sample averager given by
y(n) = Ux(n + 1) +x(n) +x(n -1)].
(3) Causal recursive filters. Implement the following difference equations. Let
NY=55.

62

INfRODUcnON TO DIGITAL SIGNAL PROCESSING

(a) y(n) = aly(n - 1) + x(n)

for a l = -0.95
(b) y(n) = aly(n - 1) - azY(n - 2) + x(n)
for a l
Let ,

= 0.95 and Wo = :n: /8.

= 2, cos( wo) and a2 = ,2.

Hint:

,n

h(n) = -.-(-) sin[(n + l)woJ

for n 2:0

SID Wo

2.5. AppHcation of a noncausal filter: a Hnear interpolator

Object. Implement a digital interpolator as a noncausal filter.
Consider the following sequence of values
0.95(n-6)

x(n)

={

for 6 :S n :S 36 and n even

otherwise

In practice, this sequence could have been obtained by placing observed data

values into the even-indexed elements of an array. The idea is to determine the
values of {x(n)} that are zeros by employing an interpolation sequence.
To perform the interpolation, let us consider the noncausal filter whose
unit-sample response is given by
h(-I) =0.5
h(O) = 1
h(l) = 0.5

Analytic results. If {y(n)} = {x(n)}

* {h(n)},

y(n) =x(n)

prove that

for n even

and

Y()
n =

x(n - 1) + x(n

+ 1)

for n odd.

Computer results. Implement the causal filter to obtain the intermediate array
{y'(n)}, then shift {y'(n)} to obtain the desired array {y(n)}. Compare the plots
of {y(n)} with {y'(n)} to verify the shift and with {x(n)} to verify the
interpolation.
Notice the strange values at the ends of the interpolated sequence. These are
artifacts produced by assuming the values of the sequence are zero outside the
specified range.
2.6. Cascading filters
Object. To implement more complicated digital filters by cascading simpler filter
sections.
Write a program to cascade NONREC and REC. To verify the filter
operation, let the input array be the delayed unit-sample sequence given by
X(N) =

{~

for N= 6,
otherwise.

For each of the filters below, perform the following:

(i) draw the filter structur.e and determine the difference equation;

DISCRETE-TIME DESCRIPTION OF SIGNALS AND SYSTEMS

63

(ii) analytically determine and sketch the unit-sample response;

(iii) compute and plot the unit-sample response;
(iv) compare the analytic and computed results_
(1) Cascade of nonrecursive filters_ Cascade
yn(n) = x(n) + 2x(n -1) + 3x(n - 2)

and

y(n) = Uyn(n

+ 1) + yn(n) + yn(n -1)]

Determine the number of nonzero elements in the unit-sample response sequence_

(2) Cascade of nonrecursive and recursive filters_ Cascade
yn(n) =x(n) - b l x(n -1)

and

for b l = r cos( wo)

yr(n) = al yr(n -1) - a2yr(n - 2) + yn(n)

Let r

= 0.95,

Wo

for a l = 2rcos(wo) and a2 =

r.

= ;r/8 and NY = 55.

(3) Cascade of recursive filters. Cascade

and
for a = 0.8 and b
2.8.)

yr(n) = ayr(n -1) +x(n)

yr'(n) = b yr'(n - 1) + yr(n)

= -0.9 and NY = 55.

(Compare this result with that of Example

1.7. DiptaI tone generators for frequency-shift keying (FSK) system.

Object. Simulate a frequency-shift keying (FSK) system by implementing two
digital tone generators at different frequencies and an appropriate switching
scheme.

In the FSK system, each bit value (lor 0) will be transformed into 32
elements of a sinusoid at a different frequency: a zero will be transmitted as
cos(won), for OSn s31, and a one will be transmitted as cos(wln), for
OSn s31.
The two frequencies will be personalized to each student: let K be the sum
modulo-8 of the last four digits of your social security number, then

wo=2;rK/32

and

WI

2;r(K + 8)/32

(if K is equal to zero, make K = 1)_

The 32 sinusoidal values are then inserted into the sequence {x(n)} to be
transmitted, the first bit sequence going into x(O) to x(31), the second valu.e into
x(32) to x(63), etc.
Analytic results. Draw a second-order recursive digital filter that acts like a digital
oscillator and generates cos(wn), for n ~ 0, for each of the two frequencies above.
Then. 32 values from the unit-sample response can then be inserted into the
transmitted sequence at the appropriate time from the desired oscillator.

Computer results
(a) Generate and plot the two 32 element sequences that correspond to the two
binary values. On the computer plot, sketch the continuous-time sinusoid from
which the samples were derived_

64

INfRODUcnON TO DIGITAL SIGNAL PROCESSING

(b) To produce an FSK sequence for future projects, write the subroutine

FSK(X)
where X (real 0lltput array) is an array of duration 96 that contains the
discrete-time sequence that results for the binary pattern 010.
Plot the output sequence. Sketch directly on the plot the continuous-time
sinusoidal function that corresponds to the samples.

FURTHER READING
Gold, Bernard, and Charles M. Rader: Digital Processing of Sig1Ulls, McGraw-Hill, 1969. [This
classic is one of the earliest books in the field.]
Oppenheim, Alan V., and Ronald W. Schafer: Digital Signal Processing, Prentice-Hall,
Englewood Cliffs, NJ, 1975. [This book was probably the most popular book in digital
signal processing, before the availability of the present text, and is suitable for advanced
study of many of the topics covered here.]
Rabiner, Lawrence R., and Bernard Gold: Theory and Application of Digital Signal Processing,
Prentice-Hall, Englewood Cliffs, NJ, 1975. [This book provides many practical aids for
digital filter design and implementation. J

REFERENCE TO TOPIC
TlDle-varying filters

Cowan, C. F. N. and P. M. Grant: Adaptive Filters, Prentice-Hall, Englewood Cliffs, NJ, 1985.
Goodwin, C. G. and K. S. Sin: Adaptive Filtering, Prediction and Control, Prentice-Hall,
Englewood Cliffs, ~J, 1984.
Haykin, Simon: Adaptive Filter Theory, Prentice-Hall, Englewood Cliffs, NJ, 1986.
Honig, M. L. and D. G. Messerschmitt: Adaptive Filters, Kluwer Academic, Boston, 1984.

CHAPTER

3
FOURIER
TRANSFORM OF
DISCRETE-TIME
SIGNALS

3.1 INTRODUcnON
In many signal processing applications, the distinguishing features of signals
and filters are most easily interpreted in the frequency domain. This chapter
describes the frequency-domain properties of discrete-time signals and the
frequency-domain behavior of discrete-time systems. The main analytic tool is
the Fourier transform of a discrete-time sequence, which can be applied either
to the unit-sample response of a filter to define the filter transfer function, or to
a data sequence to define the signal spectrum. The important properties of the
Fourier transform applied to discrete-time sequences are itemized and discussed. The magnitude and phase responses of a digital filter are examined. In
practice, discrete-time sequences are commonly obtained by sampling
continuous-time signals. The relationships between the spectra of continuoustime signals and the corresponding discrete-time sequences are investigated.
Interpolation functions to convert discrete-time sequences into continuoustime functions are also desc:ribed.

3.2 DEFINITION OF mE FOURIER

TRANSFORM
The frequency-domain analysis of signals and systems usually employs sinusoidal sequences, given by {sin(On)} and {cos(On)}, or mqre concisely {eieon },
65

66

INTRODUcnON TO DIGITAL SIGNAL PROCESSING

where w is the radian frequency. Sinusoidal sequences have the important

property that when a sinusoid of a particular frequency is applied at the input
'Of a linear time-invariant (LTI) system, then the output is also a sinusoid at the
same frequency, although it may have a different amplitude and phase. The
same statement cannot be made about other functions, such as a square-wave,
for example.
To demonstrate this form of invariance, let the.input to an LTI system be
the complex sinusoid {x(n)} = {eiaJo"}. Then, as shown in the previous chapter,
the output is equal to the convolution of the input and the unit-sample
response of the system:

y(n)=

..
..
L

h(k)x(n-k)

k=-CCI

k=--CIC

= eiaJo"

h(k) eiaJo(n-k)

h(k) e-iwok

(3.1)

k=-CID

We can interpret this last result by using complex-vector notation. Let us

define a complex-valued coefficient H(eiaJo ) that is equal to the evaluation of
the above sum, or

..

H(ei wO) =

h(k) e- iwok = IH(eiWo)1 exp{j Arg[H(eiaJo)]}

(3.2)

k--oc

Including this coefficient in Eq. (3.1), we get

y(n) = H( ei aJO) eia>o"
= IH(eiwo)1 exp{j(won

+ arg[H(eiwo)])}

(3.3)

Hence, the output sequence is a complex exponential at the same frequency

wo, but has been multiplied by a complex-valued time-invariant coefficient
H( ei aJO). This coefficient has introduced both gain and phase factors, as shown
in Fig. 3.l.
{x(nl} = {sin (won)}

wo=rr/4

{y(n)} = {iH(e l .'" llsin (won+Arg [H(e l "")])}

IH(ei"")1

Ar
012_

00

LINEAR
TIME-INVARIANT
SYS1EM

FIGURE 3-1
Magnitude and phase effects produced by a linear time-invariant system.

FOURIER TRANSFORM OF DISCRETE-TIME SIGNALS

67

This result is true for any frequency wand hence this coefficient can be
generalized to be a function of frequency. We will denote this frequency
function by H(e jW ), where
00

H(e jW ) =

h(n) e- jw"

(3.4)

This equation defines the Fourier transform of the discrete-time sequence

{h(n)}.
The Fourier transform of a sequence is said to exist if it can be expressed
in a valid functional form. Since the computation involves summing a possibly
infinite number of terms, the Fourier transform exists only for sequences that
are absolutely summable. That is, given a sequence {h(n)}, H(e jW ) exists only
when
00

(3.5)

Ih(n)l<oo

n=-CIC

When the sequence defines the unit-sample response of a digital filter,

then its Fourier transform is called the transfer function. It defines the
frequency transmission characteristics of the filter. When the sequence is a
discrete-time signal {x(n)}, then its Fourier transform is given by
00

X(e jW)

= L

x(n) e- jw"

(3.6)

n-=-oo

Then X(e jW) defines the frequency content of {x(n)}, and will be called the
signal spectrum.
Eumple 3.1. Transfer function of the 3-sample avenger. Recall that the
unit-sample response for the noncausal 3-sample averager is given by

h(n)={~ol

for-lsnsl
otherwise

Then its transfer function is given by

~

H(ei ,") =

h(n) e- J.... =

le- i ,""

n--i

= 1(ei "

+ 1 + e- i ..)

This last result can be simplified by applying the Euler identity

cia

+ e -ja = 2 COS Q'

to get

Example 3.2. Transfer function of the first-order recursive filter. Recall that the
unit-sample response is given by
h(n) = a"u(n)

for all n

68

INTRODUcnON TO DIGITAL SIGNAL PROCESSING

Because we have an infinite number of elements in {h(n)}, we must check to

determine whether it is absolutely summable:
~

Ih(n)1 =

,. _ _ 00

n-O

lain

If lal ~ 1, the sum becomes infinite and the Fourier transform does not exist. For
lal < 1, we can apply the infinite geometric sum formula

L lal"=-"=0
1-lal

for lal < 1

from which we note that the sum is finite for la I < 1. Then
~

H(eiOJ) =

L a" e-iOJn = L (a e-iOJ)"

Since la e-iOJI < 1 when lal < 1, we can apply the infinite geometric sum formula
to obtain
.
1
H(eIOJ) =
.
for lal < 1.
l-ae 1
01

The necessary condition to apply the infinite geometric sum formula in

evaluating the Fourier transform is the same for the sequence to be absolutely
summable. Hence, it is not necessary to determine whether the sequence is
absolutely summable as a separate step.
Example 3.3. Transfer
Example 2.11, we have

fuadiOD

of the secoIld-order recursive filter. From

h(n) = r" cos(woll) u(n)

Then
H(eiOJ) =

for all

11

L- r" cos(woll) e-J'"

Applying the Euler identity to the cosine, we have

i
=! [i (r e -i(OJ-"'o" + i (r e -i("'+",ol)"]
2 ,,_0
,,-0

H(eio=![i r"ei ..O"e- i.... +

r"e-iO>O"-i ....]
2 .. -0
,,-0

For Irl < 1, we can apply the infinite geometric sum formula to get

H(ei"') = U(1 -

r. e-J(W-GJol)-l + (1- r e-i(OJ+wolt1J

Combining over a common denominator, we have

.
1
2 - r(eJ....o + e-j.,O) e- i'"
H(el ....) =.
.
.
-2
...
21- r(el....O + e- IOJO) e- IOJ + r e-l~
1- r cos(I)o) e- i.,

FOURIER TRANSFORM OF DISCRETETIME SIGNALS

69

In the examples above, the sequences were written in the form of the
unit-sample response of a system {h(n)}. It is to be understood that the same
results would have been obtained for a signal sequence {x(n)}, having the
same element values. In this case, we would be calculating the spectrum,
rather than the transfer function.

3.3 IMPORTANT PROPERTIES OF THE

FOURIER TRANSFORM
We now present some important properties of the Fourier transform. These
are summarized in Table 3.1.
LINEARITY. To demonstrate that the Fourier transform is a linear transformation, we apply the superposition principle to the definition. Let
{y(n)}

= {axl(n)} + {bx2(n)}

(3.7)

Then its Fourier transform is equal to

00

Y(ei"')

= L

y(n) e- i "'"

n=-oo
00

L [axl(n) + bX2(n)] e- i "'"

n=-cc
00

=a L

00

xl(n) e- i "'" + b

n=-oo

L x2(n) e-i "'"

n=-QO

=aX1(ei ",) + bX2(e j ",)

(3.8)

Since X 1 (e J"') is the Fourier transform of {xt(n)} and X 2(ei"') is that of

{x2(n)}, we have the desired result.
PERIODICITY. The Fourier transform is a periodic function of the
continuous-valued variable w, with period 27[. This can be demonstrated by
increasing the frequency in the definition of the Fourier transform by 2k7[,
where k is any integer. We then have

H(ei ("'+2k1r

L
00

h(n) e- i (",+2kn)"

n=-ac

= L

h(n) e- i"'" e-i2knn

(3.9)

n=-tIC

For integer values of k and n, the case here, e-i2knn = 1. Hence, we have
H(ei (",+2kn = H(ei "').
This periodicity can also be noted because the variable w appears in the
function only as e i"', which is itself periodic with period 27[. If the argument of
a function is periodic, it follows that the function itself is also periodic.

70

INfRODUcnON TO DIGITAL SIGNAL PROCESSING

TABLE 3.1

Properties of the Fourier transform

1. Definition:

H(~") =

h(n) e-

(eXists when .. ~~ Ih(n)1 < 00)

I'"

2. Linearity: If

3. Periodicity:
H(~a = H(ei (OJ+2kn

for any integer k

4. Magnitude and phase functions: If

H(~OJ)

= HR(eiOJ) + jHI(eia

then
IH(~OJ)12 = Hi(~OJ) + HUe i "')
Arg[H(~OJ)

= arctan[HI(~")/ HR(eia]

5. Fourier transform of delayed sequence: For any value of no. if

(y(n)}

= {x(n -

no)}

Y(ei"') = X(eiOJ ) e-i"'''o

then

6. Fourier transform of the convolution of two sequences: If

(y(n)}

= {h(n)} * {x(n)}

then

Y(eJOJ ) = H(~"')X(ei"')

7. Fourier transform of the product of two sequences: If

{y(n)}

= {h(n)x(n)}

Y(eiOJ) =H(eiQ X(~OJ).

then

(periodic convolution)
8. For the real-valued sequence {h(n)}:
~

HR(ei"') =

h(n) cos(wn)

H1(eiOJ) = -

h(n) sin(wn)

11_-Il10

(a)
(b)
(c)
(d)
(e)
(f)
(g)

Complex-conjugate symmetry:
H(ei ") = H(ei "') = H(e- iOJ)
Real component is even function:
HR(~") = HR(e- iOJ )
Imaginary component is odd function: HI(el W ) = -H1(e- iOJ)
Magnitude spectrum is even function:
IH(~CV)I = IH(e-i"')1
Phase spectrum is odd function:
Arg[H(eiw)] = -Arg[H(e- 1W)]
If {h(n)} = {h( -n)} (even sequence), then H(~a is purely real.
If {h(n)} = {-h(-n)} (odd sequence), then H(~W) is purely imaginary.

REAL
AND
IMAGINARY
COMPONENTS
OF
THE
FOURIER
TRANSFORM. Since the Fourier transform is a complex function, it is

helpful to define alternate forms that can be plotted on a pair of twodimensional graphs. Let us consider the general case of a complex sequence
{h(n)}, which can be written as
h(n) = hR(n)

+ jhl(n)

(3.10)

FOURIER TRANSFORM OF DlSCRETETIME SIGNALS

71

where hR(n) is the real component and hJ(n) is the imaginary component. By
applying Euler's identity to the definition of the Fourier transform, we have
~

H(e i"') = ~ h(n) e-i",n = ~ [hR(n) + jhJ(n)][cos(wn) - j sin(wn)]

n=-oo

n=-~

(3.11)

This function can be expressed explicitly in terms of its real and imaginary
parts as
(3.12)

where
~

HR(ei"') = ~ [hR(n) cos(wn) + hJ(n) sin(wn)]

(3.13)

n=-oo

and
~

HJ(ei"') = ~ [hJ(n) cos(wn) - hR(n) sin(wn)]

(3.14)

n=-oo

MAGNITUDE AND
PHASE
COMPONENTS
OF THE
FOURIER
TRANSFORM. Being a complex quantity, we can also express H(eij as a

complex vector:
(3.15)
If H(ei"') denotes the transfer function of a system, it is conventional to call
IH(ei"')1 the magnitude response, and Arg[H(ei"')] the phase response. For
signal sequences, H(ei"') is the spectrum, IH(ei"')1 is called the magnitude
spectrum, and Arg[H(ei"')] is the phase spectrum. The magnitude response is
equal to
(3.16)

where H*(eiOJ) is the complex-conjugate of H(ei"'). An equivalent expression

for the magnitude is given by
IH(ei"')1 = (Hi(ei"') + Ht(ei "'ll2
(3.17)
The phase response is equal to
(3.18)

The phase is an indeterminate quantity when IH(ei"')1 =0, which occurs at

values of w for which HR(ei "') = HJ(ei"') = O.
Before giving some examples of computing the magnitude and phase
responses, we show that we can usually simplify the calculation of the Fourier
transform for real-valued sequences by considering the special properties of
their transforms.

72

INTRODUCTION TO DIGITAL SIGNAL PROCESSING

3.4 PROPERTIES OF THE FOURIER

TRANSFORM FOR REAL-VALUED
SEQUENCES
If {hen)} is a sequence of real numbers, then its Fourier transform can be
written as
00

H(e jOJ )=

00

00

L h(n)cos(wn)-j L
n=-oo

h(n)e- jOJn =

n=-tD

h(n)sin(wn)

n=-oo

(3.19)
The Fourier transform of a real sequence has the following useful symmetry
properties.
SYMMETRY OF THE REAL AND IMAGINARY PARTS. Let us consider the
effect of reversing the sign of the frequency variable in the argument of the
transform. We can then write
H(e-io = HR(e-io + jH1(e-iOJ)
(3.20)
Evaluating H(e-iO, we find
00

H(e-iO

= L

hen) e jwn =

n=-oo

00

00

L hen) cos(wn) + j L
n=-oo

hen) sin(wn)

(3.21)

n=-OC

Comparing Eqs. (3.19) and (3.21), we note that

00

HR(eiO =

hen) cos(um) = HR(e- jOJ )

(3.22)

indicating that the real part is an even function of w. Comparing the imaginary
parts, we find
00

HI(eifD) = -

L hen) sin(wn) = -H1(e-ico)

n--co

(3.23)

indicating that the imaginary part is an odd function of (J). Also, since the
complex conjugate is obtained by changing the sign of the j terms, we find
H(e- jCO ) = H*(ejO.
(3.24)
Since H(ei CO) is a periodic function of w with period 2Jr, an important
consequence of these results is that, for real-valued discrete-time sequences, all
the information in the Fourier transform is contained in the frequency range
O:s; w :s; Jr. The unspecified portion of the Fourier transform in the range
-Jr :s; w < 0 can be determined from these symmetry properties.
SYMMETRIES OF THE MAGNITUDE AND PHASE FUNCDONS. The magnitude response is also an even function of frequency, since
IH(ejOJ)1 = (H(e jOJ )H*(e jOJ
= (H*(e- jCO )H(e- jCO = IH(e-iCO)1 (3.25)

ll2

ll2

FOURIER TRANSFORM OF DISCRETETIME SIGNALS

73

The phase response can be written

Arg[H(e- iro )] = arctan[HI(e-1W)/HR(e-1W)]

= -arctan[HI(eiw)/HR(eiro)]
= -Arg[H(eiw )]

(3.26)

indicating that the phase response is an odd function of frequency. To illustrate

the above properties let us consider a few examples.
Example 3.4.
sequence. Let

Magnitude and phase spectra of the delayed unit-sample

d(n-k)={~

for n =k
otherwise

Since the unit-sample sequence is commonly used as the input to a filter, it is

considered to be a signal sequence and it is proper the call its Fourier transform a
spectrum. Let us denote the spectrum of {d(n - k)} by Dk(eIW ). Then
~

Dk(e IW ) =

2:

d(n - k) e- iwn = e- iwk

The real and imaginary parts are most easily determined by applying the Euler
identities:
and
The magnitude spectrum is equal to
IDk(eiw)1 = [cos2 (wk) + sin2 (Wk)]112 = 1
or, alternately,
IDk(eIW)1

= (e-Jwk)(eJcat)112 = 1

The phase spectrum is equal to

Arg[Dk(eiw )] = arctan[ -sin(wk)/cos(wk)] = arctan[ -tan(wk)] = -wk

The magnitude and phase spectra could also have been obtained by inspection by
noting that Dk(eJW) is already in the complex vector form. The magnitude and
phase functions are shown in Fig. 3.2 for different values of k.

In Fig. 3.2, for k = 0, the magnitude spectrum is constant with frequency

and the phase spectrum is identically zero. It is for this reason that {d(n)} is
desirable as a filter input, since the spectrum of the output sequence can then
be attributed to the filter transfer function. For the unit-sample sequence
delayed by k samples, the magnitude response is still constant with frequency,
but the phase response is linear with frequency, equal to -kw.
It is the accepted convention to plot the principal value of the phase lying
in the range [-lC, n]. Even though the phase values of +lC and -lC are
equivalent, we will allow both to occur to accommodate the plotting
conventions described below. Note from Fig. 3.2, that when the phase exceeds

74
(a)

INTRODucnON TO DIGITAL SIGNAL PROCESSING

(b)

k=O

(d")} . .

r. .

(c)

k>0

:}_.J.I_l_)

__{_d(_n-_k....

..n

-2,-1012

o1

k<O

r.

(<I(.-.))

k -10

IDu(ei"') I

-21T

-IT

If

If

I) w
2Ii

Arg[Do(ei"') J
!

-21T

-IT

t
0

-+

If

FIGURE 3-2
Magnitude and phase spectra of the unit-sample sequence. (a) Standard unit-sample sequence
(k = 0); (b) Delayed sequence (k > 0); (c) Delayed sequence (k < 0).

this range, a jump of 2.1t' is needed to bring the phase back into this range.
This principal value range is also convenient, since it is the value returned by
the commonly available ATAN2(Y, X) computer function, where Y is the
imaginary component and X is the real component.
Example 3.5. Mapitllde and pbase response of the 3-samp1e avenger. From
Example 3.1, we have

H(ei "') = !(1 + 2 cos(co

which is a real-valued function of co. The magnitude response is equal to

IH(ei"')1 = 11 +2~(CO)1
and the phase response is equal to

Arg[H(eJ"')] =

{O
1r

when H(eIOl) > 0

when H( el"') < 0

orfor - 2.1t' /3 < co < 2.1t' /3

or for -.1t'S co < -2.1t'/3
and 2.1t'/3 < co < 1r.

The magnitude and phase responses are shown in Fig. 3.3. The appropriate sign
of 1r is chosen to make the Arg[H(eio] an odd function of frequency.

As a filter, the magnitude response indicates that the 3-sample averager is

a lowpass filter. If the same discrete-time sequence were a signal, then the

FOURIER TRANSFORM OF DISCRETE-TIME SIGNALS

75

{h(n)}

...

.~"1.

-2 -10 1 2

nGURE 3-3
Magnitude and phase responses of the 3-sample averager.

magnitude spectrum indicates that it would have more energy at low

frequencies than at high frequencies.
The previous example illustrates that precautions must be taken when
determining the phase response of a filter having a real-valued transfer
function, because negative real values produce an additional phase of 1r
radians. For example, let us consider the following linear-phase form of the
transfer function that will occur repeatedly throughout the book:
H(eiQJ ) = B(eiQJ ) e- ikQJ
(3.27)
where B(eio is a real-valued function of w, called the amplitude function, that
can take on both positive and negative values. When expressing this transfer
function in the complex vector notation,
(3.28)
then the phase function must include the linear phase term and also
accommodate for the sign changes in B(ejW ). Since -1 can be expressed as
ej,., phase jumps of 1r will occur at frequencies where B(e jW ) changes sign.
This jump must be added to the phase component that is calculated from the
arctangent function. If B (eiw ) ~ 0 for all w, then no jumps will occur and the
phase function is given by - kw.
Eumple 3.6. Linear-phase Fourier transform. Let
h(O) =!

h(l) = 1

h(2) = t

h(n) = 0 otherwise

Then the Fourier transform is equal to

H(ei "') =

! + e-J'" + ~e-i2'"

76

INfRODUcnON TO DIGITAL SIGNAL PROCESSING

which can be expressed as

H(eiID) = e- iID + !e-iID(eieo + e- ieo)

= (1

+ cos(we-ieo

This last form is the linear-phase form given above. In this case, the amplitude
function is never negative and the phase response is simply - w.

PHASE JUMPS. From the previous examples, we note that there are two
occasions for which the phase function experiences discontinuities, or jumps.

1. A jump of 2Jr occurs to maintain the phase function within the principal
value range of [-n, n).
2. A jump of n occurs when B(eiev ) undergoes a change in sign.
The sign of the phase jump is chosen such that the resulting phase function is
odd and, after the jump, lies in the range [-n, n].
Example 3.7. Magnitude and phase response of the first-ordel' reausive
filter. From Example 3.2, we have H(ei"') = (1- a e- i 1 The magnitude
response is most easily obtained by calculating the squared magnitude function,
or

"'r

IH(eieoW = H(eieo)H*(e i",)

= (1- a e-i ",)-1(1- a ei ",)-l
= (1- 2a cos(w) + a 2)-1

To find the real and imaginary parts, we make the denominator real by
mUltiplying by its complex conjugate:

1
1 - a ci"'
- 1 - a e i"'1 - a ei "

H( i"') _

1- a ei "'
1- a cos(w) - ja sinew)
1 - 2a cos( w ) + a 2

Since the denominator is real, we can now separate the numerator terms to
determine the real and imaginary parts. The real part is then

/{, (ci") =
R

1- a cos(w)
1-2acos(w)+a2

The real part is an even function of frequency. The imaginary part is equal to

H1(ci"') =

-a sinew)
1- 2a cos(w) + a2

FOURIER TRANSFORM OF DISCRETETlME SIGNALS

77

a= -0.9 IH(e'''')1
a= -.0.5
I
I

---"

".

.-

/~a=-0.9

_-'---_ _ _--:.::::c:..._ _ _ _.L-.;;. w

-IT

Arg[H(e'''')]

Arg[H(e'''')]

n/2

n/2

......'"","" \

a=O.5

a=-0.5
-n/2

ftGURE 3-4
Magnitude and phase responses of first-order recursive filter for different values of Q.

and is an odd function of frequency. The phase response is then equal to

Arg[H(ei"')] = arctan[HJ(ei"')/HR(e'''')]
= arctan[ - a sine co) ]

l-acos(co)
These responses are shown in Fig. 3.4 for different values of a. When 0 < a < 1, a
lowpass filter results, and when -1 < a < 0, a highpass filter results.
Example 3.S. Magnitude and phase response of the second-order recursive
filter. From Example 33, we have
H(e i"')=

1-rcos(co~)e-I'"
1- 2r cos(coo) e- l '" + r2 e-,2O)

For simplicity of notation, let us write H(el "') in the following form
1 + a e-iO)
N(e i"')
H(el "')=
=
-l+ae- i"'+pe-,2O) D(eiw)
where N(eiw ) is the numerator and D(ei"') is the denominator. The magnitude
response is then
(N(eiO N*(e iO)lI2
.
IH(el"')1 = D(ei aJ ) D*(eIW)
_
- [1 + ~ +

[1 + 2a cos(co) + a2f12
p2 + 2a(1 + P) cps(co) + 2P coS(2CO)j112

Note from this rather complicated form that co appears only as an argument of

78

IN1"RODUcnON TO DIGrFAL SIGNAL PROCESSING

jH(eJro)j

Arg[H(elro)]
11

FIGURE 3-5
Magnitude and phase responses of second-order
recursive filter for different values of r.

the cosine function, which is an even function of w. This is reasonable, since

jH(eIW)j is itself an even function of w.
The real and imaginary parts are obtained by computing
'w

H(e

N(eiw ) D*(el "')

= D(eiw ) D*(e iW )

1 + aex + a e- iw + (af3 + ex) eiw + f3 ei2co

1 + a2 + f32 + 2ex(1 + f3) cos( w) + 2f3 cos(2w)

From this expression we find

HR(eij = 1 + aex + (a + af3 + ex) cos(w) + f3 cos(2w)

1 + a2 + f32 + 2ex(1 + f3) cos( w) + 2f3 cos(2w)
and

HI

ei )
( W

(af3

+ ex - a) sin(w) + f3 sin(2w)

= 1 + a2 + f32 + 2ex(1 + f3) cos( w) + 2f3 cos(2w)

The phase response is then given by

Ar [H(eiW)] = arctan[
(af3 + ex-a)sin(w) + f3 sin(2w)
]
g
1 + aex + (a + af3 + ex) cos(w) + f3 cos(2w)
The magnitude and phase response curves are shown in Fig. 3.5 for
various values of r.

Wo

= ,,/4 and

The magnitude response shows a sharp peak close to the frequency

which is called the resonant frequency. This peak will be
employed later to implement the passband features of a digital filter.
00

= n/4 = 000,

Example 3.9.

Magnitude and pbase response of causal3-samp1e avenger. Let

h(n)={lo

forOsns2
otherwise

FOURIER TRANSFORM OF DISCRETETIME SIGNALS

79

Then

We shall evaluate this sum using three methods: the mathematically tedious
approach, the enlightened direct approach and by applying the finite geometric
sum formula.
Mathematically tedious method. Evaluating the sum, we get
H(eiw ) = ill + e- iw + e- i2w]
= ill + cos(w) - j sin(w) + cos(2w) - j sin(2w)]

= ![l + cos(w) + cos(2w)] -1j[sin(w) + sin(2w)]

= HR(e'W) + jHI(eiw )

IH(e1W)12 = Hi(eiw) + H;(e1W)

+ cos(2w + 4 oos(w)]
= Ml + 4ooS2(W) + 400s(w)] = Ml + 2cos(wW

= ~[1 + 2(1

Taking the square root, we get

IH(e1W)1 = 111 +2cos(w)1
This result is shown in Fig. 3.6.
Direct method. To allow an Euler identity substitution, we factor out e- 1W :
H(eiw) = !e-iw(eiw + 1 + e- iw )
=

W + 2cos(we- iw = B(eiw)e-iw

From which we get

IH(e1W)1 = ill + 200s(w)1

h(n)

~ III

-.....l~O~1-72-!j-4~-

....
~ n

Arg(H(eJO>)]

:~

-1T

FlGURE 3-6
Magnitude and phase responses of causal 3-sample averager.

80

INTRODUCfION TO DIGITAL SIGNAL PROCESSING

and

Arg[H(eJW )] =

{-W
-Wlr

for -2lr/3<w<2lr/3,
for -lr:S W < -2lr/3

or when B(eJW ) >0

or when B(eJW ) < O.

and 2lr/3 < w:S lr,

The phase is undefined at frequency points where IH(eJW)1 = o.
Geometric sum method. Applying the finite geometric sum formula with
r = e- Jw, we get

The above form can be simplified by substituting the following identity:

1 - e-'" = e- j<>i2( e J<>i2 - e- Ja/2)
= 2j e- J<>i2 sin(a-/2)

After this substitution, we get

e - j3 0.>'2(e j3 0.>'2 _ e - J3 0.>'2)

H(e JW ) = .J

3 e- 0.>'2(e J0.>'2 - e J0.>'2)

1

= 3

e -jw sin(3w/2)
sin(w/2)

B( j"') -jw
e e

This last form is the amplitude with linear-phase form described above in Eq.
(3.27). From that discussion, we find
IH(eJW)1 = Isin(3w/2)1
13 sin(w/2)1
and

Arg[H(eIW )] =

{-W

-w lr

when B(eIW ) > 0

when B(e jW ) <0

The magnitude response, an even function of w, is expressed as the ratio of

two odd functions of w. The equivalence of this form with that obtained with the
previous two methods is left as an exercise.

3.5 PROGRAM TO EVALUATE THE

FOURIER TRANSFORM BY COMPUTER
To provide computer verification of the analysis above, we will describe a
simple program to evaluate the Fourier transform of a real-valued sequence at
N p points over the frequency range 0:::; W < 2n. In the next chapter, we will
present a more efficient but mathematically more complicated algorithm for
performing this procedure.
Although not the most efficient method, the evaluation of the real and
imaginary parts at a frequency point Wk can easily be performed directly from
the definition. For example, let h(n), for 0:::; n:::; NH - 1, be a sequence of

FOURIER TRANSFORM OF DlSCRETETIME SIGNALS

81

real-valued elements. Then, the real part of the Fourier transform evaluated at
frequency Wk> denoted by HR(e iWk ), is given by
MI-l

HR(eJWk )

2:

(3.29)

hen) cos(wkn)

n=O

(3.30)

To obtain a discrete-frequency approximation to the continuous-frequency

function H(e iw ), we can evaluate the real and imaginary parts at a set of N p
equally-spaced points over 0:::; W < 2n, i.e. at Wk = 2kn 1Np , for 0:::; k :::;
Np-l.
When processing arrays whose indexing starts at one, rather than zero,
the above equations for the real and imaginary parts must be modified. If
H(K), for 1:::; K :::;NH, is the array of discrete-time values, with H(1) equal to
h(O), then the array containing the real part, denot"Ai by HR(K), is given by
HR(K) =

NH

2:

H(N) * cos(2n(K - 1)(N - 1)1 N p )

for1:::;K:::;Np

(3.31)

N=l

The frequency of the Kth cosine sequence is equal to 2n(K -1)INp , for
1:::; K:::; N p The array containing the imaginary part, HI(K), is
HI(K) = -

NH

2:

H(N) *sin(2n(K --1)(N -l)INp )

for 1 :::;K:::;Np

(3.32)

N=l

Using the programming style described in the pre lious chapter, the
computation of the Fourier transform should be performed with the subprogram, denoted by

FT(H,NH,HR,HI,NP)
where
H

(real input array) contains the elements of the discrete-time

sequence;
NH
(integer input) is the number of elements in H;
HR
(real output array) contains the real part of the Fourier transfonn;
HI
(real output array) contains the imaginary part of the Fourier
transform;
NP
(integer input) is the number of frequency points over O:s W < 2n.
For displaying the Fourier transform, the magnitude and phase responses
are conventionally used. The magnitude and phase responses should be
computed using the subprogram denoted by

MAGPHS(HR,HI,HMAG,HPHS,NP)

82

INlRODUcnON TO DIGITAL SIGNAL PROCESSING

where
HR
HI
HMAG
HPHS
NP

(input real array) contains the real part of the Fourier

transform;
(input real array) contains the imaginary part of the Fourier
transform;
(output real array) contains the magnitude values~
(output real array) contains the phase values; and
(input integer) is the number of points over O:s OJ < 2Jr.

The array of magnitude values can be computed as

HMAG(K) = SQRT(HR(K) * HR(K) + HI(K) * HI(K)]

for 1 :s K:s NP
(3.33)

The array of phase values can be computed as

HPHS(K) = ATAN2(HI(K), HR(K

for 1 :s K :s NP

(3.34)

When both HI(K) and HR(K) are zero, the arctangent function should not be
used, but rather a value of zero should be assigned for the pha~e value.
Some programming languages only have the ATAN(Y IX) function,
which generates the phase in the first and fourth quadrants, i.e. in the range
[-n/2, Jr/2]. It is left as an exercise for the student to determine the principal
value of the phase using the ATAN function and the sign of the real part.
Because of the symmetry properties for real-valued sequences, we need
to evaluate the plot the Fourier transform only over the range 0 :s OJ :S Jr. Since
both the frequency points OJ = 0 and OJ = Jr are included, the number of points
to be computed is equal to one more than half the points covering one entire
period, or (Np /2) + 1.

3.6 USE OF THE FOURIER

TRANSFORM IN SIGNAL PROCESSING
The popularity of the Fourier transform for digital signal processing is due to
the fact that the complicated convolutional operation in the time domain i~
converted to a much simpler multiplicative operation in the frequency domain.
To demonstrate this, let {yen)} be the output of a LTI system. Then, its
Fourier transform is given by
Y(eio

yen) e-i "'"

(3.35)

Applying the convolutional definition of the output, we have

Y(eiO

= ,,~'" L~

h(k)x(n - k)] e- i"'"

(3.36)

FOURIER TRANSFORM OF DISCRETE-TIME SIGNALS

MUltiplying this last equation by unity, in the form eJwk e- Jwk , we get
'"
.,
Y(eiw ) =
[h(k) e-iwkx(n - k) e-Jw(n-k)]

2: 2:

83

(3.37)

n=-oo k=-oo

Since h(k) e- iwk is not a function of n, it can be factored out of the summation
with respect to n, giving
Y(eiw )

= k~'" [h(k) e- iwk n~'" x(n -

k) e-iW(n-k)]

(3.38)

In the sum over n, the term (n - k) represents a shift of {x(n)} by k samples.

Further, note that the value of the exponent follows the index of {x(n)}. Since
the summation is over the entire sequence, this shift is inconsequential in
evaluating the sum. By replacing (n - k) with the dummy variable m, we note
that the second summation becomes the definition of X(e iw ), the Fourier
transform of {x(n)}. Since X(e iaJ ) is not a function of k, we can factor it out of
the summation over k. The remaining terms in the sum then define the Fourier
transform of {hen)}. Hence, we have arrived at the desired result
(3.39)

This last result indicates that the input spectrum X(e iW ) is changed through a
multiplicative operation by the filter transfer function H(e iW ) to produce the
output spectrum Y(ei aJ ). Applying the complex vector notation to the output
spectrum, we obtain
Y(eiaJ) = IY(eiaJ)1 exp{j Arg{Y(ei"')]}
= IH(ei"')1 exp{j Arg[H(ei"')]} IX(eiaJ)1 exp{j Arg[X(ei"')]}

(3.40)

=IH(ei"')IIX(ei"')1 exp{j(Arg[H(eiaJ)] + Arg[X(ei"')])}

Equating the magnitudes and the phases, we find
(3.41)

and
(3.42)
In words, the output magnitude spectrum is equal to the product of the input
magnitude spectrum and the filter magnitude response. The output phase
spectrum is equal to the sum of the input phase spectrum and the filter phase
response.
Example 3.10. Sipal processing with the Fourier transform. Let the input
sequence be given by
x(n) =

{~

for -1 Sn S 1

otherwise

84 INTRODUCTION TO DIGITAL SIGNAL PROCESSING

and the filter unit-sample response by

a"

h(n) = { 0

for n 2:0
otherwise

From Examples 3.1 and 3.2, we have

X(e iw ) = i{1 + 2cos(w
and

The output spectrum is then equal to

Y(eiw) =H(eiw)X(eiW ) =

~(; ::o;~~

The output squared-magnitude spectrum is equal to

IY(eiw)12 = (1 + 2 cOS(~)(1 + 2 COS(W)
3(1 - a e IW) 3(1 - a elW )

_
(1+2cos(wW
- 9(1- 2a cos(w) + a2 )
The output phase spectrum can be computed in two ways. First, it can be
obtained from the real and imaginary components of Y(ei W ) , determined from
iW

Y(eiw ) = (1 + 2 COS(~)(1- a e )
3(1- a e- IW ) 1- a elw
_ (1- a cos w)(1 + 2 cos(w - ja sin(w)(1 + 2cos(w
3(1- 2a cos(w) + a2 )
Then
j ..
_
(-aSin(w)(1+2COS(W)_
(-asin(W)
Arg[Y(e )]-arctan (1-acos(w(1+2cos(w -arctan 1-acos(w)

The phase can also be determined by adding the phase response of the filter to
the phase spectrum of the input:
Arg[Y(e i "')] = Arg[H(eIW)] + Arg[X(e i "')]

-a sin(w) )

= 0 + arctan ( 1- a cos(w)

3.7 FOURIER TRANSFORM OF SPECIAL

SEQUENCES
In this section, we consider the Fourier transform of several special sequences
that simplify the computations. These results will be helpful for the analysis of
digital filters later.

FOURIER TRANSFORM OF DlSCRETETIME SIGNALS

8S

FOURIER TRANSFORM OF ADELA YED SEQUENCE. Let us examine the

effect on the Fourier transform when a sequence is delayed by k sample points,
as shown in Fig. 3.7. Let us denote the delayed sequence by {hk{n)} =
{h(n - k)}. Then its Fourier transform, denoted by Hk{e jCJJ ), is equal to
Hk(e jCJJ )

2:'"

hk{n) e- jOJn

n=-oo

'"
2:
h{n n=-OCI

k) e- jWft

(3.43)

(h(n)}

Original
sequence

... llil I
o

Delayed
sequence

"..

1 2

(hJn)}

........
II!II.;.
o
1 k-l k k+l

Arg[H(e JW )]

IH(el"')1

1T

-It
.-

CI)

I.-

I
I

It

CI)

-.-

-------1

I
I

F1GURE

~7

Comparison of the magnitude and phase spectra for original and delayed sequences.

CI)

86

INTRODUcnON TO DIGITAL SIGNAL PROCESSING

Multiplying the right side by e- iwk e1Wk, we have

Hk(e jW ) =

00

2:

hen - k) e-1W(n-k) e- jwk

n=-IXI

=e- jwk

00

hen - k) e-iw(n-k)

(3.44)

n=-oo

By substituting m = n - k in the above summation, we recognize the definition

of H(e iw ). Hence,
Hk(e iw ) = H(eiw) e- iwk
(3.45)
By expressing this result in complex vector form, we note
IHk(eiw)1

= IH(eiW)1

(3.46)

and
(3.47)

In words, when we delay a sequence by k samples, where k is any positive or

negative integer, the magnitude of the Fourier transform is unchanged. The
phase is altered by adding a linear-with-frequency phase term whose slope is
equal to -k. The utility of this result will become evident in computing the
Fourier transform of symmetric sequences in the next section.
FOURIER TRANSFORM OF SEQUENCES WITH SYMMETRIES. In the design

of digital filters, we commonly encounter sequences that have a symmetry

about some time point, denoted by Ns . We consider sequences that have even
symmetry and those that have odd symmetry separately.
Even-symmetric sequences are those whose elements about the point Ns
have equal value, or
heNs + m) = heNs - m)

for all m

(3.48)

The two cases, depending on whether Ns is an integer, are shown in Fig. 3.8.
The special case when Ns = 0 defines an even sequence, for which
hen) = h( -n)

for all n

(3.49)

Antisymmetric sequences are those whose elements about the point Ns are
of equal magnitude, but of opposite sign, or
heNs + n) = -heNs - n)

for all n

(3.50)

The two cases, depending on whether Ns is an integer, are shown in Fig. 3.8.
When Ns is an integer, heNs) must equal zero for {hen)} to be an
antisymmetric sequence. The special case when Ns = 0 defines an odd
sequence, for which
hen) = -h(-n)

for all n

(3.51)

FOURIER TRANSFORM OF DISCRETE-TIME SIGNALS

ANTISYMMETRIC SEQUENCES

EVEN-SYMMETRIC SEQUENCES

11 INs TIT

NsINlEGER

Ns NONINlEGER

87

~J-LT-",---o.>

-.-J'LT1-...1-1-1,1c'-1

Ns

nGURE 3-8
Symmetric and antisymmetric sequences shown for integer and noninteger values of symmetry
point N s -

Fourier transform of even-symmetric sequences. Let us first consider the

Fourier transform of an even sequence, followed by that of a general
even-symmetric sequence. Starting with the definition of the Fourier transform, we have
H(eiw ) =

..

h(n) e- iWrl

n=-oo
-1

..

h(n) e- iwn

n=-oo

+ h(O) + L h(n) e- jWrl

n=1

Letting m = -n in the first sum, we obtain

H(e IOJ )

..

=L

h( -m) el<mt

m=1

h(n) e- iwn

n=1

Noting that h(n) = h( -n), we have

H(e iOJ )

..

+ h(O) + L

..

= h(O) + L

h(n)(e- iwn

+ elOJn )

n=1

= h(O)

+2 L

h(n) cos(wn)

(3.54)

n=1

Hence, the Fourier transform of an even sequence is a real-valued function of

frequency. The magnitude function is then
.
{H(eIW )
when H(eiw ) ~O
(3.55)
IH(eIOJ)1 = -H(eiw )
when H(eiOJ ) <0

88

INTRODUcnON TO DlOITAL SIGNAL PROCESSING

Since its imaginary part is identically zero, the phase function computed from
the arctangent function is zero. The phase response then only accommodates
the sign change of H(e iw ):

Arg[H(ei W )] = {O
Jr,

for
for

iw

(J)
(J)

over which H(e ) > 0

over which H(e iW ) < 0

(3.56)

Recall that the phase is not defined when H(e iW) = O.

Example 3.11. Fourier transform of an even sequence. Let
h(n) = 1

for -2:5 n :52,

=0 otherwise

Then its Fourier transform is given by

n--2

= 1 + 2 cos( w)

+ 2 cos(2w)

Any even-symmetric sequence can be obtained from an even sequence by

applying a suitable delay. The Fourier transform of. a delayed sequence was
shown above to have the same magnitude function as the original sequence,
but also to contain an additionallinear-with-frequency component in the phase
function. To illustrate this in a more general fashion, let us consider the even
sequence {hE(n)}, which is defined for -MS;ns;M, as shown in Fig. 3.9. If
{hE(n)} is the unit-sample response of a digital filter, then the filter is
noncausal because of the nonzero elements for n < o. If a causal filter having
the same magnitude response is desired, we need only add sufficient delay so
that the first nonzero element occurs at n = O. As shown in Fig. 3.9, we can
obtain a causal symmetric sequence {hs(n)} by delaying {hE(n)} by M
samples, or
(3.57)
for OS;n s;2M

{~(II)}
. .I

I
I III_
IIII

-M

I . .

-lOI2 M

{1Is(n)} = {h.!(n-M)}

..I

I II

I I I I II I .

-1 01

2M

FIGURE 3-9
n

Oclaying an even sequence to produce a causal symmetric sequen<:e.

FOURIER TRANSFORM OF DISCRETE-TIME SIGNALS

89

Then, if Hs(eJeo)is the Fourier transform of {hs(n)} and HE(eJW ) is that of

{hE(e jW )}, we have
(3.58)
and
(3.59)
Example 3.U. Fourier transform of an even-symmetric sequence. Let
h(n) = 1

for 0:5n:54

This is just the delayed version of the sequence given in Example 3.11. We can
get the Fourier transform by simply adding the linear phase function - 2w to the
result of Example 3.11 to include the delay.
We verify this result by directly computing the Fourier transform of the
above sequence, which is given by
4

H(e jW ) =

L e-

...

-0

Applying the finite geometric sum formula, we get

1- e- Jsw

H(e JUJ ) = - -JUJ

l_e- -

e - jS on ejSwi2 _ e - jSwi2
e ion

eiwi2 _

e iwi2

= sin(5w/2) e- J2UJ

sin(w/2)
The phase is a linear function of frequency with slope equal to -2. It is left as an
exercise for the student to show that the two amplitude expressions are
equivalent.

Recall that an even sequence has zero phase, except for the phase jumps
of .;rr caused by the sign reversal of the amplitude of the Fourier transform. If
these jumps are ignored, we are left with a phase that is linear with frequency,
which will be called a linear-phase function. This linear-phase property will be
important for filter design in the following chapters.

FOURIER TRANSFORM OF ANflSYMMETRIC SEQUENCES. We first consider the Fourier transform of an odd sequence, followed by that for the
general antisymmetric sequence. For an odd sequence, h(O) = 0, and the
Fourier transform can be written as
00

H(e iw ) = ~ h(n)[e-iwn - eieM ]

n=l
00

=-

j2 ~ h(n) sin(wn)
n=l

(3.60)

90

INTRODUCTION TO DIGITAL SIGNAL PROCESSING

Hence, the Fourier transform of an odd sequence of n is an imaginary function

of frequency. Since only sine functions appear in the sum, H(e iw ) is an odd
function of 0, i.e. H(ei W ) = -H(e- 1CO ). Since the real part is identically zero,
and j = ei J<i2, we can write the function in terms of its imaginary part as

H(e iCO ) = jHI(eiCO ) = HI(eico) ei J<i2

The magnitude response is then equal to

IH(eiCO)1

= HI(e iCO )

for HI(eico ) > 0

=-H1(eico )

for HI(eico ) < 0

(3.61)

and the phase is equal to

for
for

Ar [H(eico) = {lC/2
g
lC/2 lC

0
0

over which HI(ei CO) > 0

over which HI(eico) < 0

(3.62)

Example 3.13. Maguitude and phase response of aD odd sequence. Let

-I
h(n) =

01

forn=-2
for n =2
othe~

The Fourier transform is equal to

H(e 1tu) = _e j2tu + e- j2tu = -2j sin(2w) = 2sin(2w) e-j"fl
The imaginary part is then equal to
H,( e itu) = -2 sin(2w)
The magnitude spectrum is given by
IH(e'"')1 = 12 sin(2w)1
and the phase spectrum is
Arg[H(el."')] = {lC/2

-lC/2

for H,(ei tu ) > 0,

for H,( eitu) < 0,

or for lC/2 < w < lC

or for 0 < w < lC/2

The magnitude and phase functions are shown in Fig. 3.10.

For odd sequences, a jump of lC degrees will always occur at 0 = O. This

happens because HI(eico) is an odd function of 0, which means that a sign
reversal in the amplitude function must occur at 0 = O.
An antisymmetric sequence can be obtained from an odd sequence by
applying a delay. Let us consider the odd sequence {ho(n)}, which is defined
for -M::::: n::::: M. We can obtain a causal antisymmetric sequence {hA(n)} by
delaying {ho(n)} by M samples:
hA(n) = ho(n - M)

for O:::::n :::::2M

(3.63)

Then, we have
(3.64)

fOURIER TRANSFORM OF DISCRETE-TIME SIGNALS

91

IH(eI'")1

'VV\
o

1tfl

1t

-,.

Arg[H(eJ")
1t.
1tfll-

I
1t

-lTfl
-IT'-

FIGURE 3-10
Magnitude and phase spectra for an odd sequence.

and
(3.65)

Ignoring the phase jumps due to sign reversals in the amplitude function,
including the one at (J) = 0, we are also left with a linear-phase function for
antisymmetric sequences. This result is true even when M is not an integer,
although a strictly odd sequence does not exist in this case. This is illustrated in
the next example.
Example 3.14. Fourier transform of aD aDtisymmetric sequence. In this example, we will cOnsider the case when Ns is not an integer. As shown in Fig. 3.11,
IH(e)")1
{h(n)}

-2 -1

~,

-IT

1t

FIGURE 3-11
Magnitude and phase spectra
for an antisymmetric sequence.

92

INTRODUCfION TO DIGITAL SIGNAL PROCESSING

let h(O) = -!, h(1) = 1, h(2) = -1, h(3) =! and h(n) = 0, otherwise. Then, the
Fourier transform is given by
H(ei"') = _! + e-J'" - e-J2", + !e- i3 ",
= -!(1- e- i3 ",) + (e- j'" - e- j2 ",)
= _e-i3WI2H(ei30i2 _ e- i30i2) + ei0i2 _ e- i0i2]
= -je- i3",I2(sin(3ro/2) + 2 sin(ro/2

Substituting -j = e- iJ1f2 , we finally get

H(e iW ) = (sin(3rol2) + 2 sin(ro/2e- i(3wu)/2

When Ns is not an integer, the phase spectrum still has a linear slope equal to
-Nsro. In terms of Eq. (3.65), M = ~.

3.8 THE INVERSE FOURIER

TRANSFORM
The inverse Fourier transform allows us to determine the unit-sample response
of a filter for which we know only the transfer function. Similarly, it also helps
us determine the output discrete-time sequence of a filter from the output
spectrum.
The Fourier transform H(e iw) is a periodic function of the continuousvalued variable (i), with period 2.1l'. Hence, it can be expanded as a Fourier
series. We show that the coefficients of the series terms are equal to the values
of the discrete-time sequence {h(n)}. We first state the Fourier series

argument:
Let s(t) be a periodic function of the continuous-valued variable t with
period T. Then s(t) can be expressed in the following Fourier series
expansion:
co

s(t) =

an e- iot

(3.66)

where On = 2.1l'n/T is the nth harmonic frequency, and the coefficients

are given by

an

=-1 fTI2
T

-TI2

s(t) eiO.,t dt

(3.67)

In exactly the same fashion, H(eiw) can be expressed as the following Fourier
series expansion
co

H(eiw )
where 8n

= 2.1l'n/2.1l' = n,

L en e- iBw
n=-co

(3.68)

can be thought of as a harmonic frequency,.and

(3.69)

FOURIER TRANSFORM OF DISCRETE-TIME SIGNALS

93

Note the similarity of the definition of the Fourier transform in Eq. (3.4) and
the Fourier series expansion of H(eiw ) given in Eq. (3.68). This similarity
allows us to define the inverse Fourier transform as
(3.70)
This integral solution for the inverse Fourier transform is useful for analytic
purposes, as shown later, but it is usually very difficult to evaluate for typical
functional forms of H(eiw ). An alternate, and occasionally more helpful,
method of determining the values of h(n) follows directly from the definition
of the Fourier transform. equal to
00

H(eiw ) =

2:

h(n) e-IW>l

n=-oo

If H(eJW) can be expressed as a series of complex exponentials, then h(n) is

simply the coefficient of e- jOft This is illustrated in the following example.
EDIIlple 3.ts. Inyerse Fourier transform of the 3-sample ayenger. Let
H(e'''') =

W + 2cos(w

Employing the Euler identity for the cosine we obtain

H(ei"') = i(1 + eiw + e- i"').
Noting that h(n) is the coefficient of e- i... , we have by inspection
h(O)=h(I)=h(-I)=t
h(n) = 0
forn<-landn>1
Alternately, we could perform the integration. The value of h(n) is then
h(n) = - 1
2.n=

ilf 3(1 + e'''' + e-'''') e' .... dw

1

_If

J.. ("

[ei'" + eiw(II+I) + e1w(,,-I)] dw

6.n-L.
_J.. (el -

- 6.nsin(n.n-)
3n.n-

e- i.... e i1i(n+l) - e-i 1i(" +I) e i1i(n-l) - e- i1i(n-I)

jn
+
j(n + I)
+
j(n -1)

= --+

sin[(n + 1).n-J sin[(n -1).n-J

+ ---='-'---~
3(n + 1).n3(n -1).n-

Recall that sin(n.n-) =0 for integer values of n. When the denominator is also
zero, an indeterminate form results. Applying L'Hospital's argument, we find
that

lim sin(E)/E = 1
.-0

94

INTRODUcnON TO DIGITAL SIGNAL PROCESSING

Applying this approach to each of the three terms above, we obtain

hen) =

{t

for n = -1, 0,1

otherwise

H H(e iw ) is a rational function of e-i ""', then a series expansion can be

obtained either by performing a division or by computing the Taylor series

expansion, as illustrated in Example 3.16.
Example 3.16. Inverse Fourier transform for the first-order recursive filter. Let
H(e jOJ) = (1- a e-JOJ)-l
Computing the Taylor series expansion, we have
!:f(ejOJ) = 1 + a e- jOJ + a 2 e- j2 Q) + a 3 e- j3OJ + ... + a k e- Jwk + ...
Equating the coefficient of e- jOlk with h(k), we have

a"

h(n)= { 0,

for n 2:0
otherwise

Example 3.17 illustrates the application of frequency domain analysis to

signal processing, in which the output sequence of a filter is determined from
the output spectrum.
Example 3.17. DetermiDing the output sequence from
spectrum. The output spectrum of Example 3.10 is given by
Y(ei OJ)

the

output

=! ejQ) + 1 + ~-j..,
3 1-ae-JOJ

Expanding each term in the parenthesis in a Taylor series expansion, we get

Y(ejO/) = l[ejO/ + a + a 2 e-jO/ + a 3 e- j2 0/ + a 4 e- j3w + .. .
+ 1 + a e-jQ) + a 2 e- j2 0/ + a 3 e- j3 Q) + .. .
+

e-JO/+a e-j2OJ+a2e-J30/+ ... ]

Combining all the terms having the same complex exponential factor and
equating h(k) with the coefficient of e -jOlk, the output sequence is found to be
yen) = l[an+1u(n + 1) + a"u(n) + an-1u(n - 1)]

3.9 FOURIER TRANSFORM OF THE

PRODUCT OF TWO DISCRETE-TIME
SEQUENCES
In the design of nonrecursive filters, we commonly encounter the product of
two sequences {w(n)} and {h(n)}. Let the sequence of term-by-term products

FOURIER TRANSFORM OF DISCRETETIME SIGNALS

95

be denoted by {v(n)}, where

v(n) = w(n) h(n)

for all n

(3.71)

w(n) h(n) e- icm

(3.72)

Then its Fourier transform is given by

V(ei lil) =

"=-00
To determine a more useful form for V(eiW ), we first express h(n) as an
inverse Fourier transform given by
(3.73)
We thus get
(3.74)
Interchanging the order of summation and integration, we have

(i

V(eiO=~

(w(n)e-i(W-cr)n)H(eiJda
2Jf )-_ 11_-"

(3.75)

The term in the parentheses is similar to the definition of the Fourier transform
of {w(n)}, but at a frequency argument of co - a, or W(ei("'-cr. Including this
Fourier transform into the integral, we get
V(eio =

~ (- H(ei J W(ei(o>-cr da
2Jf 1-_

= H(eiO

W(ei"')

(3.76)

This form is that of a convolution operation involving two functions of the

con~uous-valued variable co. Because both these functions are periodic, the
integral is evaluated over one period. We denoted this operation as a circular
convolution. Hence, we have shown that multiplication in the time domain
corresponds to convolution in the .frequency domain. We will use ibis result in
the next section to determine the Fourier transform of sampled continuoustime data.
Eu.pIe 3.18. Fourier traDsfolDl of the prodac:t of two sequaaces. The most
common application of this result occurs in truncating infinite-duration sequences,
as illustrated here.
Let us consider the infinite duration sequence given by
h(n) -

Sin!:;~4)

for all II,

The corresponding H(eI"') is equal to the ideallowpass filter characteristic shown

Fig. 3.12. To truncate this sequence to -3 s n s 3, we multiply h(n) by the

in

96

INTRODUcnON TO DIGrrAL SIGNAL PROCESSING

... 011. 'II

OJ

rlrlhl. 11' .'''....

1~]lllllL. .
(~" I ,,")hi"))

II r

.. .I-3 I 0

-3

HC'")

It. ..
3

FIGURE 3-12
Multiplication of two discrete-time sequences with the corresponding spectral convolution
operation.

following sequence, called a window function,

w(n) =

{~

for -3::s n ::s 3

otherwise

The windowed sequence is then

v(n)

= w(n) h(n) =sin(:rn/4)/(:rn/4)

for -3::s n::s 3

V(ei "'),

we must first compute W(e;") and

and is shown in Fig. 3.12. To find
perform the convolution. The Fourier transform of {wen)} is equal to
3

W(d"') =

2:

e- j .... = ePO>

n--3

. 1 - e- j7 .,
=eI3 ..
.
1 - e I'"

2: e-;'"
n-O

sin(3w/2)
sin (w /2)

The result of the convolution is shown in Fig. 3.12. The main effect of the
convGlution is that the sharp transitions in H(ej .,) have been made smoother. We
will explore this multiplication with a window sequence in more detail in the
design of finite-impulse response filters in a later chapter.

3.10 SAMPLING A CONfINUOUS

FUNcrION TO GENERATE A
SEQUENCE
In practice, discrete-time sequences are often obtained by sampling
continuous-time signals. Typically, this occurs when signals are to be processed

FOURIER TRANSFORM OF DISCRETETIME SIGNALS

97

by computer; situations include speech analysis and extracting information

from sensors for bioengineering and robotic applications. We will apply these
results to sampling the Fourier transform in the next chapter. An application
of sampling continuous-time signals will also occur in a later chapter in the
design of digital filters. In this section we relate the spectrum of the
continuous-time signal to that of the discrete-time sequence that is obtained by
sampling.
We proceed by first stating the relationship between the continuous-time
and discrete-time spectra and then provide the rather involved mathematical
proof. Let us define the continuous-time, or analog, signal by CA(t). Its
spectrum is then given by CAUQ), where Q is the continuous-time frequency,
having units radians per second. For continuous-time signals, the Fourier
transform is defined by
(3.77)

which is similar to the definition of the Fourier transform of discrete-time

signals, but with the sum replaced by an integral. When CA(t) is sampled every
T. seconds, we obtain the discrete-time sequence {cA(nT.)}, whose values are
equal to
(3.78)

T. is commonly called the sampling period and has units seconds per sample
interval. The discrete-time Fourier transform for such sampled data sequences
is defined by
(3.79)
n=-oo

where w has units radians per sample interval. This discrete-time spectrum is
related to the continuous-time spectrum CA(jQ) by
C(e

iw

= ~ k~OO CAOw/T. + j2kn/T.)

(3.80)

In words, the spectrum of the discrete-time sequence, derived by sampling a

continuous-time function, is equal to the sum of an infinite set of shifted
continuous-time spectra. A graphical interpretation of this identity is shown in
Fig. 3.13.
Proof. To show this relationship, we take the following approach. First, we
treat sampling as a multiplicative operation involving the continuous-time
signal CA(t) and a train of Dirac delta functions occurring every T. seconds.

This multiplication produces the intermediate function of continuous-time

cs(t). The continuous-time spectrum of cs(t), denoted by CsOQ), is shown to
have the same analytic form as Cs(eiW ). Finally, we determine CsOQ) in terms
of the original analog spectrum CA(jQ).

98

INfRODUcnON TO DIGITAL SIGNAL PROCESSING

f\
-:bE
Ts

-It

Ts

7t

2r.

Ts

Ts

co,n

FIGURE 3-13
Spectrum of a discrete-time sequence obtained by sampling a continuous-time signal.

We first define the Dirac delta function and present its properties that are
important for this application. By definition, the Dirac delta function, denoted
by c5(t), has infinite height and zero width, but its area, or weight, is equal to
unity. To illustrate how we can approach a Dirac delta function, consider the
rectangular function r.(t), whose width is E and whose height is lIE, so its area
is equal to one, as shown in Fig. 3.14. We can define c5(t) as
c5(t) = lim rE(t)

(3.81)

....0

The following three properties of the Dirac delta function are useful
here.
Property 1. The product of a function and a delta function. In multiplying any continuous-time function f(t), the weight of the delta function
becomes equal to the value of the function at the instant of the delta function,
or
f(t) c5(t - to)

= f(t o) c5(t -

to)

(3.82)

Property 2. The integral of the product (The sifting property). The

integral of the product of any function f(t) and a delta function is a constant,
equal to the value of the function at the instant of the delta function. To show
this, we integrate the product to give
fJ(t) c5(t - to) dt = f./(t o) c5(t - to) dt

,.(1)

II

-'h

'h

o
1"==0

FIGURE 3-14
Limiting proCedure to obtain a Dirac delta function.

(3.83)

FOURIER TRANSFORM OF DISCRETETIME SIGNALS

99

by Property 1. Since f(t o) is not a function of t, we can take it out of the

integral to give
f(t o) L"'", o(t - to) dt = f(t o).

This last equality follows because the area of o(t) is unity. The Dirac delta
function is used to select, or sift, the value of the function at t = to
Property 3. Convolving a function with a delta function (The shifting
property). When the function f(t) is convolved with o(t - to), the origin of
f(t) is shift to t = to, or
o(t - to) * f(t)

= f(t -

to)

To show this, we write out the continuous-time convolution equation,

(3.84)

['" 0(.\ - to)f(t -.\) cIA

and define the variable =.\ - t. We then have

[",f(-) 0(<< + t - to) d = f(t - to)

(3.85)

This last result follows by Property 2.

An infinite train of Dirac delta functions occurring every
denoted by or;(t) and can be expressed as
or,(t) =

'"

L o(t n=-oo

nT.)

T. seconds is
(3.86)

We can now use these properties of the Dirac delta function to determine
the relationship between the spectra. The sampling process is modeled by the
multiplication operation given by
for all

= cA(t)

L'"

o(t - nT.)

for all t

(3.87)

n=-CXl

as shown in Fig. 3.15. Since cs(t) is a continuous-time function, its spectrum

Cs(jQ) is given by
. Cs(jQ) = ['" cs(t) e-

[J

cA(t)

jQt

dt

n~'" o(t -

(3.88)

nT.)] e-

jQt

dt

(3.89)

Changing the order of integration and summation, and applying the sifting

100

INTRODUcnON TO DIGITAL SIGNAL PROCESSING

FIGURE 3-15
Sampling a continuous-time signal can be treated as a multiplicative process with a train of Dirac
impulses.

property of the Dirac delta function, we have

n~.., [f.., 6(t -

Cs(jQ) =

nT.) CA(t) e-

jDt

dt]

...

2:

cA(nT.) e-jnCT,

(3.90)

We have arrived at our first intermediate result. By

(3.79), we observe

compari~g

Eq. (3.90) with

(3.91)

That is, the form of C(eiQ) can be obtained by determining CsOQ) and
replacing Q by (J) I T.. The division of (J) by T. scales the values of (J) to be the
same as Q. For sampled-data applications, it is conventional to express w in
the range [-niTs, nlT.), to have the same numerical values as Q. The units of
wiT. are radians per second, just as for the continuous-time frequency.
To complete the proof, we must determine Cs(jQ) in terms of the
continuous-time spectrum CA(jQ). Since 6r.(t) is a periodic function, it can be
expressed in the following Fourier series expansion given by
..,

6 T ,(t) =
where Q n

= 23m1T., and

i
T.

2: an e-jc
n=-oo

T
,12

an = -1

(3.92)

6 T,(t) ejQ.r dt

(3.93)

-T,12

Over the interval - T./2 :$ t < T./2, the impulse train 6r.(t) is represented by
only the single Dirac impulse situated at t = 0, or 6(t). Making this substitution
and applying the sifting property, we get

i
T.

1
an =-

T
,12

-T,12

1 j o=1
6(t) dQ.rdt=-e
.

T.

T.

for all n

(3.94)

FOURIER TRANSFORM OF DISCRETE-TIME SIGNALS

101

Substituting this Fourier series expansion into the sampling process, we get

cs(t) = CA(t) . ar,(t)

1
= CA(t) . -

=-

oo
L

L
00

T. n=-oo

Tan=-cc

e- IQ t

CA(t) e-IQ t

(3.95) .

The continuous-time spectrum of cs(t) is then equal to

CsGQ) = roo Cs(t) e-

jQt

dt
(3.%)

Interchanging the order of summation and integration, we obtain

Cs(jQ) =

~ n~oo [f"oo CA(t) e-j(Q+Q)t dt ]

(3.97)

Evaluating the kth term of the sum, denoted by CS,k(jQ), we find

CS,k(jQ) =

L:

CA(t) e-j(Q+Qk)t dt

= CAOQ + jQk)

(3.98)

This result is simply the original analog spectrum that has been shifted from
0=0 to 0 = -Ok, as shown in Fig. 3.16. The total spectrum is then the sum
of these shifted versions and is given by
(3.99)

(= -2ITk)
T,

FIGURE 3-16
Shifted version of the analog spectrum.

102

INTRODucnON TO DIGITAL SIGNAL PROCESSING

Cs(jQ) is a periodic function of Q, with period 21f/T.. More concisely, Cs(jQ)

is called the periodic extension of CAOQ).
.
This last identity is the second desired intermediate result. Combining
this result with that given by (3.91), we arrive at the final desired identity given
in (3.80).
(Q.E.D.)
NYQUIST CRITERION. When sampling a continuous-time signal cA(t) to
produce the sequence {cA(nT.)}, we want to ensure that all the information in
the original signal is retained in the samples. There will be no information loss
if we can exactly recover the continuous-time signal from the samples. To
determine the condition under which there is no information loss, let us
consider CA(t) to have a bandlimited spectrum, or one for which
(3.100)
as shown in Fig. 3.17(a). As demonstrated above, when CA(t) is sampled with
sampling period T., then the spectrum of the sampled signal CsOQ) is the
periodic extension of CA(jQ) with period 21f/T., as shown in Fig. 3.17(b). The

(b)

T, < 7T/OM

C (j11)

Ideallowpass filter

_"' ' _~. .J.:- ~-,- r- I<: ~ZS:_- _-" "*-_-_- -:-! -~-L-i-LL_~'_'-+~
- 211

IS

-11
--

Ts

11"

11"

-11

Ts

211

Ts

FIGURE 317
Relationship between continuous-time spectrum of a signal and the spectrum of the discrete-time
sequence obtained by sampling the signal with sampling period T.. (a) Original spectrum of
continuous-time signal; (b) spectrum of sampled sequence when 11M < nIT.; (c) spectrum of
sampled sequence when Q M > nIT.. The latter case illustrates the aliasing error.

FOURIER TRANSFORM OF DISCRETE-TIME SIGNALS

-m

-1[

Ts

IS

1[

21[

IS

IS

CAoy,

-1[

1[

IS

IS

103

~(jn)

U
-

IS

- n
Ts

FIGURE 3-18
Two candidates for the continuoustime spectrum when aliasing occurs_

form of Cs(jQ) in the frequency range [-.nIT., .nIT] is identical to CA(jQ) if

or

(3.101)

In this case, there is no overlap in the adjacent spectral components and we

can then extract CA(jQ) exactly from Cs(jQ) by employing the ideal lowpass
filter magnitude response shown in Fig. 3.17(b).
If T. is chosen to be greater than .n/QM , spectral overlap occurs in the
periodic extension, given by Eq. (3.99), and the form of Cs(jQ) in the range
-.nIT. $. Q < .nIT. is then no longer similar to CAOQ) as shown in Fig. 3.17(c).
This overlap, caused by sampling at too Iowa rate, produces an irretrievable
error in the spectral values, called aliasing. The true spectral shape is
irretrievable since many different CAOQ) functions can produce the same
Cs(jQ). Two possible candidates are shown in Fig. 3.18.
A useful sampling measure is the sampling rate, Is = liT.. Recalling that
Q = 2.n1, and defining 1M as the highest frequency component in the signal,
then no spectral overlap will occur if
(3.102)
To prevent aliasing error, more than two samples are required per period of
the highest frequency sinusoidal component present in the signal. The smallest
sampling rate before aliasing occurs for a particular continuous-time signal is
called the Nyquist rate. The effect of aliasing is most easily demonstrated by
considering a sinusoidal signal in Example 3.19.
Example 3.19. Demonstration of aliasing. Let cA(t) = cos(Qot), for all t, with
Qo = lOOn. This function is sampled to produce the sequence
cs(n 7;) = cos(nQo7;)

for all n

104

INTRODUCl10N TO DIGITAL SIGNAL PROCESSING

(b)

(a)

CA(j!l)

o 'lo

It

T,

k' J
\

-It

T,

~21T

/ 'lo

2!..
T,

C,(eJOl)

C,(el"')

-21T

I
0

2.. - 'loT 0 'loT IT

~~\----~--~~;IT

-+'loT
T,

I T~IT I

21T

--'loT - +'loT
T,
T,

--------

FIGURE 3-19
Demonstration of aliasing caused by choosing too large. (a) Unaliased version, when < 11:/0.0 ;
(b) aliasing caused when > 11:/0.0 In this latter case a low-frequency cosine function is apparent
when the sequence is interpolated.

r.

r.

r.

The continuous-time spectrum of CA(t) is equal to two Dirac delta functions

as shown in Fig. 3.19. Sampling CA(t) causes the spectrum to be extended

periodically:
Cs(jO) =

2;',,~~ [6(0 -

0 0 - 2nK/T.) + 6(0 + 0 0 - 2nKIT.)]

To demonstrate aliasing we note the position of the impulses situated

within -KIT. ~ 0 < KIT., as we vary the sampling period T.. The sampling period
dictated by the sampling rate is T. < KIOo(=10-2s).
For T. < KIOo, say T. = 10- 3 s, the delta functions falling within -KIT. ~
0< KIT., or ....:l000K ~ 0 < 1000K, are only those for the n = 0 term above.
These occur at 0 = Oo. In the discrete-time domain, these impulses will
produce the desired sequence {cos(nOoT.)}.
For T.> KIOo, say T. = 1.1 X 10- 2 s, the impulse in the range 0< 0 < KIT.,
or 0 < 0 < 9OK, is at
0= (2KIT.) - 0 0 = 180K -lOOK = 80K rad S-1

from the term 6(0 + 0 0 - 2KIT.), or the k = 1 term. The impulse located in the
range -KIT. < 0 ~ 0 is at
0= 0 0 - 1K/T. = lOOK -180n = -BOn rad S-1

FOURIER TRANSFORM OF DISCRETE-TIME SIGNALS

105

from b(g - go + 2nllT.), or the k = -1 term_ In this case, the sampled sequence
appears to be derived from a continuous-time function at a lower frequency
2n/T. - go = BOn rad S-I_ This is the aliasing phenomenon_

This effect of aliasing errors is occasionally observed in films of moving

cars in which the wheels appear to be turning in the direction opposite to that
expected from the motion of the car.
ANfI-ALIASING FILTER. In practice, the frequency range contammg sig-

nificant power in the original continuous-time signal may be larger than that
containing the desired information. This commonly occurs when a lowfrequency signal is corrupted by high-frequency noise or in modulated signals,
such as broadcast radio. If sampled at the rate dictated by the Nyquist criterion
for the desired analog signal, unwanted high frequency signals would cause
aliasing errors to occur.
To prevent aliasing errors caused by these undesired high-frequency
signals, an analog lowpass filter, called an anti-aliasing filter, must be used.
This filter is applied to the analog signal prior to sampling and reduces the
power in the analog signal for the frequency range beyond Q = niT.. In
practice, the spectral magnitude level for Q ~ JrIT. should be less than 1%
(-40 dB) of the important spectral features in the desired signal spectrum to
prevent significant aliasing.
For example, let the observed analog signal contain power in an
arbitrarily large frequency range, but have relevant information only in the
frequency range -Q R $ Q $ QR. We wish to determine the appropriate
sampling period T.. The Nyquist criterion tells us that T. must be less than
Jr/QR. If sampled at this rate, the higher frequency components in the analog
signal will be aliased into the relevant discrete-time signal. The anti-aliasing
filter must satisfy two conditions:
Condition 1. The analog signal components with frequencies less than
must be negligibly attenuated by the filter.
Condition 2. To prevent aliasing in the relevant analog frequency range,
the components for IQI ~2JrIT. - QR must be attenuated. After sampling, the
components in this frequency range fall into the range -Q R $ Q $ Q R when
the periodic extension is formed.
QR

We consider how these conditions can be met in Chapter 8, where the design
of analog lowpass filters is discussed.

3.11 RECONSTRUcnON OF
CONTINUOUSTIME SIGNALS FROM
DISCRETETIME SEQUENCES
In many practical applications, discrete-time sequences are required to be
transformed into continuous-time signals. One example is the digital audio

106

lNTRODUcnON TO DIGITAL SIGNAL PROCESSING

CS<jO)

D
27t

-T,

/:\ D
I

_~-OM

T,

OM

7t

27t

T,

T,

H,(jO)
I 1""1- - - . ; . ; . : . . . . . ; . . . . . - - - ,
1
I

L
0

I
I
I

7t

o
FIGURE 3-20
The spectrum of an analog signal can be obtained by multiplying the spectrum of the sample values
with the transfer function of an ideal analog lowpass filter. In practice, H,(jQ) must be
approximated.

compact disk system, in which the audio signal is stored as digital samples on
the disk and must ultimately be converted into the analog signal that drive the
speakers. Another is the digital controller in a robot system, that uses the
digital computer to process the sensor data, but requires an analog signal to
control the actuator motors. In these cases, we are faced with the converse
problem from that of sampling CA(t) to obtain {cA(nT.)}. The relevant
question now becomes how to recover the analog signal CA(t) from the
sampled sequence {csCnT.)}.
We begin by considering the unaliased spectrum of Cs(jQ), shown in Fig.
3.20. CA(jQ) then has the same form as Cs(jQ) over -rclT.:5 Q < rclT.. Let us
consider an analog lowpass filter with the following ideal frequency response
for -rclT.:5Q<rclT.
otherwise

(3.103)

We could retrieve CA(jQ) from Cs(jQ) by filtering with an ideallowpass filter,

or
for all Q

(3.104)

To interpret this result in the time domain, let us examine the corresponding
convolution equation given by

CA(t) =

L"'", cs(c)hICt -

c) dc

(3.105)

The function csCt) is the C<?ntinuous-time sampled signal composed of a series

<

FOURIER TRANSFORM OF DISCRETETIME SIGNALS

107

Approximation to C,(t)

"/

...

FIGURE 3-21
Approximating the sequence of weighted Dirac
delta functions with a sequence of narrow pulses
having different amplitudes.

_ _~--l~~_--':'-_u....._. n T,

of weighted Dirac delta functions. In practice, this continuous-time function

could be approximated by narrow pulses whose amplitudes are proportional to
the corresponding sample values cs(nT.), as shown in Fig. 3.21- The function
hI(t) is the continuous-time impulse response of the ideallowpass filter, found
by the inverse Fourier transform given by

=----j2m

~ Sinc(tlT.)

1 sinemIT.)
lttlT.

= T.

(3.106)

The Sinc(x) function is shown in Fig. 3.22. The function is equal to zero for
integral values of the sampling time, i.e. t = nT., where n is a non-zero-valued
integer. For t = 0, we find
hI(O) = lim.!. sin(mlT.)
1-+0

T.

:Jrt I T.

= .!.lim (mIT.) T. 1-+0

(1/3!)(lttlT.)3 + ...
Jrt I T.

lIT.

(3.107)

Because cs(t) is composed of a series of weighted Dirac delta functions,

the convolution produces a sum of shifted and scaled versions of hI(t), as
demonstrated by Property 3 of the Dirac delta function. As shown in Fig. 3.23,

~.
-3T-2

-T

2T 3T

4T

FIGURE 3-22
The Sinc(tIT) function, used
for interpolation, is the impulse response of the ideal
analog lowpass filter.

188

INTRODUcnON TO DIGITAL SIGNAL PROCESSING

ftGURE 3-23
In the time domain, the interpolation
process is treated as the convolution of
weighted and delayed Sinc(t/T) functions.
The values of the sequence {cs(nT.)} are
shown as dots. The continuous-time
function is shown in heavy solid line.

the original analog signal can be reconstructed from the samples by

oo

CA(t) =

L cs(nT.) hI(t n=-co

=.!.

T.n=-oo

nT.)

cs(nT.) sin(.n-[t - nT.]IT.)

.n-[t-nT.]IT.

(3.108)

The action of hI(t) in this convolution is to fill in, or interpolate, the

values of the continuous-time function between the sample points. Such
functions are called interpolation functions. Because of the infinite limits of the
Sinc(x) function, it is awkward to apply practically. In practice, the ideal
interpolating filter must be approximated by a causal lowpass analog filter.
PRACI1CAL INfERPOLATION FILTERS. The problem with the ideal interpolation function hI(t) is that it extends in time from minus infinity to plus
infinity. Hence, it incorporates the entire sequence, all the future as well as all
the past samples, into the computation of the interpolated function at any time
point. Especially in control situations, we cannot usually wait to observe the
entire sequence before an interpolation is performed. In this section we
consider some of the properties that practical interpolators should have and
provide examples of two commonly employed interpolators.
Let us consider the impulse response of a practical continuous-time
interpolating filter to be given by </>(t). In practice </>(t) should have the
following three properties:
(1)

(2)

(3)

</>(0) = 1
</>(nT.) = 0

for n *0

foo I</>(t)1 dt <

(3.109)

00

The first two properties permit the interpolated function to equal the sample
values at the sampling instants. Property (3) guarantees that the interpolated
function remains finite between the sampling instants. Example 3.20 considers
two commonly employed interpolation functions that have these properties.

FOURIER TRANSFORM OF DlSCRETETIME SIGNALS

109

1;

~(t)

-T,

nGURE 3-24

1;

Two commonly used interpolation functions.

Eumple 3.20. Zero-order hold and linear iDterpolaton. The impulse response
of the zero-order hold interpolator 4>o(t) is given by

4>o(t) =

{~

forO~t<T.

otherwise

That of the linear interpolator 4>1(t) is given by

for -T~t~T
otherwise
These interpolating functions are shown in Fig. 3.24. Notice that both these
functions have the three desirable properties.
Examples of interpolated functions are shown in Fig. 3.25. The zero-order
hold is very easy to implement by quickly charging and discharging a capacitor at
the sampling instants, but the interpolated function CAO(t) exhibits discontinuities

T.

nGURE 3-25
Results of applying the two interpolation functions to the same
sequence.

110

INTRODUcnON TO DIGITAL SIGNAL PROCESSING

at the sampling times. These discontinuities may be filtered by the mechanical

inertia of the syst!,!m that is driven with the signal.
The linear interpolator produces a smoother function CA1 (t). But this
interpolation is noncausal, since it must look ahead to determine the slope value.
To make it causal, a delay can be introduced so that the interpolated function lags
behind the samples by one sampling period. In control systems, this delay is
crucial and may restrict the use of the linear interpolator.

3.12 SUMMARY
In this chapter, we presented the frequency-domain procedures for analyzing
digital filters. The Fourier transform of the unit-sample response provided the
frequency filtering characteristics of the digital filter, while the Fourier
transform of a data sequence indicated its frequency content. The Fourier
transform was shown to be a periodic function of w, with period 2Jr, and for
real-valued sequences to have a magnitude function that is an even function of
frequency and a phase function that is odd. This allows the entire function to
be determined from the values for O:s; w :s; Jr. Even or odd real-valued
sequences were shown to have purely real or imaginary Fourier transforms,
while sequences with symmetries have linear phase functions. The operations
of convolution and multiplication were shown to be converted into the other
when the Fourier transform was computed. The inverse Fourier transform,was
defined through a Fourier series argument, but found to be complicated to
apply directly. An alternate method involving a Taylor expansion was found to
be simpler, but did not provide an analytic form for the time series. The
relationship between continuous-time and discrete-time frequency functions
was discussed and so were procedures to convert a continuous-time function
into a discrete-time sequence. This chapter ended with a short discussion on
reconstructing the analog signal from the sampled sequence by using an
interpolating filter. In the next chapter, we will investigate the signal
processing implications of evaluating the Fourier transform at a finite number
of points by computer.

FURmER READING
Antoniou, Andreas: Digital Filters: Analysis and Design, McGraw-Hill, New York, 1979.
Bracewell, R. N.: The Fourier Transform and Its Applications, 2d ed., McGraw-Hill, New York,
1978.
Childers, D. G. (ed.): Modern Spectral Analysis, IEEE Press, New York, 1978. [This is a book of
reprints which presents a historical perspective of the evolution of spectral analysis.]
Gold, Bernard and Charles M. Rader: Digital Processing of Signals, McGraw-Hill, New York,
1969.
Hamming, R. W.: Digital Filters, 2d ed., Prentice-Hall, Englewood Cliffs, NJ, 1983.
Jong, M. T.: Methods of Discrete Signal and System Analysis, McGraw-Hill, New York, 1982.
Oetkin, G., T. W. Parks, and H. W. Schussler: "New results in the design of digital
interpolators", IEEE Trans. ASSP ASSp23: 301-309 (1975).
Oppenheim, Alan, V., and Ronald W. Schafer: Digital Signal Processing, Prentice-Hall,
Englewood Cliffs, NJ, 1975.

FOURIER TRANSFORM OF DISCRETETIME SIGNALS

111

Papoulis, A.: The Fourier Integral and its Applications, McGraw-Hill, New York, 1962.
Papoulis, A.: Signal Analysis, McGraw-Hili, New York, 1977. [This book contains the
continuous-time results referred to in this chapter.)
Peled, A. and B. Liu: Digital Signal Processing: Theory, Design and Implementation. Wiley, New
York,1976.
PhiIlips, C. L. and'H. T, Nagle, Jr.: Digital Control Systems Analysis and Design, Prentice-Hall,
Englewood Cliffs, NJ, 1984.
Rabiner, Lawrence R. and Bernard Gold: Theory and Application of Digital Signal Processing,
Prentice-Hall, Englewood Cliffs, NJ, 1975.
Schwartz, Mischa; Information, Transmission, Modulation and Noise, McGraw-Hill, New York,
1980. [This book contains the continuous-time results referred to in this chapter.)
Schwartz, Mischa, and Leonard Shaw: Signal Processing: Discrete Spectral Analysis, Detection and
Estimation, McGraw-Hill, New York, 1975.

PROBLEMS
3.1. Given the following unit-sample response sequences compute the Fourier
transform and sketch the magnitude and phase functions.
(a) h( -1) = 1, h(l) = 1, h(n) =0, otherwise.
(b) h(O) = 1, h(l) = -1, hen) = 0, otherwise.
(c) h( -1) = 1, h(O) = -2, h(l) = 1, hen) = 0, otherwise.
3.2. Determine whether the following sequences have Fourier transforms and
compute the Fourier transform when it exists.
(a) hen) = 1, for n ~O, = 0, otherwise.
(b) h(n) = 1, 0:5 n :5 N, = 0, otherwise.
(c) hen) = 1, for all n.
(d) h(n)=1, for -N:5n:5N, =0, otherwise.
3.3. An important theorem of signal processing is Parseva/'s theorem, that can be
stated as
~
1
"~~ x(n)x*(n)=2n -n X(e1W)X*(e1W)dw

fn

Prove this theorem by substituting the definition of X(e 1W ) into the integral and
evaluating.
3.4. Show that a causal filter that has more than two nonzero elements in its
unit-sample response cannot have a phase response that is zero for all
frequencies.
3.S. Consider the sequence {g(n)} that is the time-reversed version of {hen)}, or
g(n)=h(-n)

for all n

(a) Let the H(e1W ) be the transfer function of the filter whose unit-sample
response is {hen)}. Compute the transfer function G(eiW ) in terms of H(e 1W ).
(b) Consider taking the output of the filter described by H(eiw ) and letting it be
the input to a second filter having the transfer function G(e iw). The total
transfer function is then G(e iW) H(e~'). What are the magnitude and phase
responses of this cascade combination?
3.6. For the following difference equations, determine the transfer function of the
filter and sketch the hlagnitude and phase responses.
(a) yen) = x(n) - x(n - N). Let N = 4.
(b) yen) = 2r cos(O) yen - 1) + r2y(n - 2) + x(n). Let 0 = n/4 and r = 0.9.

ill

INTRODUcnON TO DIGITAL SIGNAL PROCESSING

3.7. For the following magnitude and phase functions, determine and sketch the
unit-sample response, the digital filter implementation and the difference
equation.
(a) IH(ei"')1 = [20 + 16 cos(2w WI2,
j

Arg[H(e ,")] = arctan (2

:S:~~;2 )

(b) IH(ej'")1 = [17 - 8 coS(2W)]II2,

4 sin(2w) )
Arg[H(el."')] = arctan ( --~~
1-4cos(2w)
(c) IH(ej"')1 = [1.81-1.8cos(w)]1I2,
. ]
( 0.9sin(w) )
Arg[H (e1"') = -w + arctan 1 _ 0.9 cos(2w)
3.8. Show that the two different analytic forms for the magnitude response in Example
3.9 are equivalent.

3.9. Some computer languages provide only the ATAN(Y/X) function, which
produces the radian angles in the range - : /2 :5 w :5 : /2. Describe a procedure
to determine the principal value of the arctangent over the range -::5 w < :,
from the result of the ATAN(Y/X) function and the values of X and Y.
3.10. If we wanted to generate the sequence
x(n) = cos(won/lO),

for all n

by sampling the continuous-time signal

X(I) = COS(Ool)

for all 1

where Wo = 0 0 1;. What is the value of the sampling period 1;?

3.11. Exploring tbe aliasing phenomenon. The follQwing continuous-time Fourier
transform pairs for the sinusoids of infinite duration are given by:
if

U(I) = COS(Ool)

then

U(jO) = (:12)[c5(0 - 0 0 ) + c5(Q'+ 0 0 )]

if

V(I) = sin(Ool)

then

V(jO) = (j:12)[ c5(0 - 0 0 )

and
-

c5(0 + 0 0 )]

Determine and sketch the real and imaginary parts. of the spectra Us(jO) and
Vs(jO) for the continuous-time sampled functions for the following three sampling
periods:
(i) T. < : 100
(ii) 1;>>>:100
(iii) 1; = :100
3.U. SampUng. bandlimited signal. Let CA(I) have a spectrum CA(jO) that is nonzero
only over the frequency range 0 1 :5101:5 O 2 , as shown in Fig. P3.12. It can be
shown that cA(I) can be sampled at a rate that has to be greater than twice the
bandwidth, O 2 - 0 1 , rather than twice the highest frequency O2 , and yet have no
information loss. This can be demonstrated by doing the following steps.
(a) Let 1; = 2:1(0 2 - 0 1), Determine and sketch the periodically extended
spectrum of the time series.

FOURIER TRANSFORM OF DISCRETETIME SIGNALS

f'1

---.. .uo,....----'o,....---:!!--~U-~---.;.~n
-fi2 -fi]

fi]

fi2

113

FIGURE P3-U
Bandpass spectrum of continuoustime signal.

(b) Sketch the magnitude response of the bandpass filter that retrieves the
desired signal spectrum.
(c) The impulse response response of this filter acts as the interpolation function.
Assume the magnitude response of part (b) and a zero phase response to find
the impulse response. Sketch this continuous-time function. Consider
whether this waveform is reasonable by considering its values at the sampling
instants.

COMPUTER PROJECTS
3.1. High frequency sequences
Objed. Observe the shape of sinusoidal sequences for frequencies in the range
:re/2 < (J) <:re, since they can often appear to be confusing.
Computer results. Use TESTO to obtain plots of either the sine or cosine sequence
for the following values of frequency: (J) = :re/2, 3:re/4, 7:re/8.
Aoalytic results. Determine and sketch the corresponding continuous-time sinusoid that connects the dots of the computer plot.
3.2. Computer evaluation of the Fourier transform
Objed. Write and verify the program that evaluates the Fourier transform of a
real-valued sequence at a specified number of frequency points in the range
Os (J) < 21&.
Analytic results. Determine and sketch the magnitude and phase spectra of the
sequences given below and compare them with the computer results. In all cases,
evaluate the sum in the Fourier transform.
Computer results. Write the subprogram

FT(H.NH.HR.HI.NFT)
where H is the input array containing the real-valued discrete-time sequence of
length NH, HR and HI are the output arrays of size NFT containing the real and
imaginary values of the Fourier transform, with the Kth element corresponding to
the frequency point (J)K = 21&(K - 1)/NFT, for 1 S K S NFT. For this project, let
NFT = 128, which will allow comparison with the fast Fourier transform in the
next chapter.
A clever way of obtaining the value of :re in the program is by using the
arctangent function: PI = 4. *ATAN(1.).
To compute the magnitude and phase spectral sequences, write

MAGPHS(HR. HI. HMAG. HPHS.NFT)

114

INTRODUCTION TO DIGITAL SIGNAL PROCESSING

which converts the input real HR and imaginary HI arrays of the Ff of size NFf
into the output magnitude HMAG array and phase HPHS array. Use the
FORTRAN function ATAN2(XI, XR), where XI is the imaginary component
and XR is the real component, to calculate the principal value of the phase, i.e.,
in the range [-:re, :re). Undefined values of the phase, when HR(K) = HI(K) = 0,
should be set to zero.
To plot the magnitude and phase sequences, use the PLOTMR and
PLOTPR subroutines given in Appendix A. These plot the 65 elements of
HMAG and HPHS that correspond to the range 0!5 W !5 :re.
Compute the spectra for the following sequences.
(a) h(O) = t h(l) =~, h(2) = 1, h(3) = t h(4) = land h(n) = 0, otherwise.
(b) The unit-sample response of first-order recursive filter for the following
values of a: 0.8, 0.95, -0.8, and -0.95 (compute first 128 elements).
3.3. Becoming familiar with the phase spectrum
Object. Illustrate the jumps in the phase response curves.
Analytic results. Determine and sketch the phase spectra of the sequences given
below and compare them with the computer results. In all cases, evaluate the sum
in the Fourier transform. Indicate the phase jump location and type on "the
sketches and computer plots.
Computer project. Use the programs Ff and MAGPHS written in Project 3.2 to
compute the phase and plot with the subroutine PLOTPR given in Appendix A.
Use NFf = 128.
(a) Jumps of 2:rc in the phase spectrum. Compute the phase spectrum of the
delayed unit-sample sequence {d(n - k)} for k = 2. Repeat for k = 3 and
k = 10. Sketch the linear phase functions on the plots.
(b) Jumps of :re in the phase spectrum. Compute the phase spectrum of the
following sequence: h(O) = h(l) = h(2) = 1, and h(n) = 0, otherwise. /
(c) Both jumps of :re and 2:re. Compute the phase spectrum for the following
sequence: h(n) = 1 for 1 !5 n !5 5, and h(n) = 0, otherwise. Indicate the causes
'Of the jumps.
3.4. Demonstration of aliasing
Object. Demonstrate the aliasing phenomenon both analytically and numerically.
Analytic results. Do Problem 3.11 and compare the results with those produced
by the computer.

Computer results. For this project use the TESTO.FOR program in Appendix A
to generate both the sine and cosine sequences. Rather than varying the sampling
period with respect to the frequency of the sinusoid, we will keep the sampling
period constant, equal to 1, and vary the frequency. Demonstrate the aliasing
effect by plotting the time-domain sequences which result when the frequency 0 0
is varied to correspond to the three cases in Problem 3.11.
3.S. Frequency domain interpretation of interpolating sequences
Object. Explore the frequency-domain behavior of two discrete-time interpolation
functions.
"
Consider the sequence
v(n) = cos( won)

for 0!5 n s;NV -1

FOURIER TRANSFORM OF DISCRETE-TIME SIGNALS

115

V(e JW )

---O~--~W~o--------------------~n--+W

X(e JW )

FIGURE Proj 3-5

o
with

Cl.Io=

IT

Spectra for Project 3.5_

n/16 and NY = 128 and x(n) such that

x(n) =v(n)
x(n) =0

for n even, 0:5 n :5 NY - 1

for n odd, O:5n:5NY - L

The idea is to determine the values of {x(n)} that were zeroed, by employing
interpolation functions as was done in ProJect 2.5_
The spectra of the two sequences, denoted by V(ei"') and X(eiOJ ) are shown
in Fig_ Proj. 3.5. They are related by
~

X(eIOJ ) =

2:

V(ei(OJ+kn-

k=-oo

In words, X(eiOJ) is the periodic extension of V(eiOJ), but with period n.

AaaIytk results. We will compare tlte performance of two discrete-time
interpolators:
Zero-order interpolator, where

for0:5n:51
otherwise
First-order interpolator, where

0.5
1.0
ht(n) =

~.5

forn=-1
forn=O
forn=1
otherwise

116

INTRODucnON TO DIGITAL SIGNAL PROCESSING

(a) Determine the magnitude responses IHo(eiw)1 and IH1(e jW )I. Sketch both
responses on the same graph for 0:5 W :5 n.
(b) Let {yo(n)} and {Yl(n)} be the outputs of the zero-order and first-order
interpolators, respectively, when {x(n)} is the input. Compute the magnitude
spectra IYo(eiw)1 and 1y'(e'W)I. Sketch both spectra on the same graph.
Computer results. The interpolation cap be implemented as a nonrecursive filter:

SUBROUTINE NONREC (X, NX, BCOEF, NB, Y)

where X is the input array of size NX, BCOEF is the coefficient array of size NB,
and Y is the output array of size NX + NCOEF - 1. For the nonrecursive filter
recall that the coefficient sequence is equal to the unit-sample response. Since
NX = 128, dimension the Y array sufficiently large to contain the output
sequence.
(a) Generate {v(n)} and {x(n)} for wo=n/16. Plot the first 55 samples.
Implement the two interpolators and plot 55 samples of the output sequences
{yo(n)} and {Yl(n)}. The first-order interpolator corresponds to a noncausal
filter. Repeat the above for Wo = n / 4.
(b) Compute and plot the magnitude spectra of {v(n)} and {x(n)}. Use
NFf= 128. Do the same for Yo(n) and Yl(n). for O:5n :5127.
(c) Explain why there is a difference in the interpolator performance for the two
values of woo
3.6. Detection fiben for the FSK system
Object. Implement and investigate the filters that differentiate between the two
sequences generated by the FSK subroutine written in project in Chapter 2.
In the FSK subroutine, 32 samples of two different frequencies were used to code
the binary values 0 and 1 to be transmitted over communication lines. In this
project, the samples of the two sinusoids will serve a second purpose, other than
the code. They will correspond to the unit-sample response of the two receiver
filters, but in time-reversed order. This is, if {ho(n)} is the unit-sample response of
the "zero-frequency (Ko) filter". then
ho(n) = cos[2nKo(31 - n )/32]

for O:5n :531

The unit-sample response of the "one-frequency (K 1 ) filter" is given by

for O:5n :531
When the unit-sample response is matched in this manner to the desired input
signal, the filter is called a matched filter.
Analytic re:mIts
(a) Sketch the block diagram of the interpolation filter as a nonrecursive
structure, i.e. containing only feed-forward coefficients.
(b) Determine the magnitude response of each filter. (Hint: use the finite
geometric sum formula.)

Computer results
(a) Compute and plot the magnitude spectra of the two 32-sample sequences (use
NFf = 128). Use PLOTMR to plot only the first 65 points. Indicate the
location of the two FSK frequency values on each magnitude spectrum plot.

FOURIER TRANSFORM OF DISCRETETIME SIGNALS

117

(b) Implement the two filters using the NONREC subprogram. Compute and
plot the magnitude response of the two filters (use NFl' = 128). Compare
these to part a.
(c) Generate and plot the four discrete-time sequences that are obtained by
passing each of the 32-bit sequences into each of the two filters.
(d) Generate and plot the discrete-time output sequence for each filter when the
input is the FSK sequence corresponding to bit pattern 010. Based on these
output patterns, describe a simple procedure to determine which binary value
was transmitted.

CHAPTER

.4
THE DISCRETE
FOURIER
TRANSFORM

4.1 INTRODUCl10N
The discrete Fourier transform (DFr) sequence, denoted by {H(k)}, allows us
to evaluate the Fourier transform H(e i "') on a digital computer. This complexvalued sequence is obtained by sampling one period of the Fourier transform
at a finite number of frequency points. The DFf is important for two reasons.
First, it allows us to determine the frequency content of a signal, that is, to
perform spectral analysis. The second application of the DFr is to perform
filtering operations in the frequency domain. After defining the DFr and the
inverse discrete Fourier transform (IDFf), we examine some of its properties
that are relevant for signal processing. The DFr can be applied directly to
causal, finite-duration sequences. Difficulties arise, however, when we try to
apply it to noncausal or infinite-duration sequences. These difficulties are
examined and solutions are suggested. The. fast Fourier transform (FFr)
algorithm is presented as an efficient computational procedure for evaluating
the DFr. Its efficiency is the main reason much current signal processing is
performed in the frequency domain. The FFr algorithm is described in both an
intuitive and analytical form using the decimation-in-time approach.

us

119

THE DISCRETE FOURIER TRANSFORM

4.2 THE DEFINITION OF THE

DISCRETE FOURIER TRANSFORM (DFT)
The DFf is the finite-duration discrete-frequency sequence that is obtained by
sampling one period of the Fourier transform. This sampling is conventionally
performed at N equally spaced points over the period extending over
Os w <2K, or at
forOsksN-1

(4.1)

If {h(n)} is a discrete-time sequence with the Fourier transform H(eiw), then

the DFf, denoted by {H(k)}, is defined as

forOsksN-1

(4.2)

The DFf sequence starts at k = 0, corresponding to w = 0, but does not

include k = N, corresponding to w = 2K.
Let us consider the consequences of sampling H(eirD). In the previous
chapter, it was shown that the Fourier transform H(eij is periodic in w, with
period 2K. Its 'Fourier series expansion was then computed and the coefficients
were found to be equal to the discrete-time sequence {h(n)}. It was also
shown that when a function of continuous-valued time is sampled with
sampling period To, then the spectrum of the resulting discrete-time sequence
becomes a periodic function of frequency with period 2K/T.. A similar result
can be shown to be true here. As shown in Fig. 4.1, when H(ei"'). is sampled
with sampling period w. = 2K/N, the corresponding discrete-time sequence,
denoted by {Ji(n)}, becomes periodic in time with period 2K/w. = N. This

H(el"')

h(n)

/~.'

......... 111 .........

-1 01

21r(J)

tbysamPling

H(k)

h{n)

.111 ..... 111 ..... 111

...
N
-N

III 1 , , I III.

-1

0 1 2

N-IN

FIGURE 4-1
A periodic extension of a time sequence results when a spectrum is sampled. Case for N
samples of spectrum over 0 S w < 2.n-. Bracket indicates Npoint DFr sequence.

=8

120

INfRODUcnON TO DIGITAL SIGNAL PROCESSING

periodic discrete-time sequence can be expressed in terms of {hen)} as

co

hen) =

hen + mN)

(4.3)

m=-oo

The sequence {hen)} is called the periodic extension of {hen)}. The number of
sample points in one period of the spectrum N is also equal to the period of
{hen)}.
The importance of this result is apparent when we employ the DFf for
signal processing and filtering in the frequency domain, for it is {hen)} that is
the corresponding sequence that is being used in the time domain. This chapter
describes the implications and pitfalls that accompany this result.
We are free to choose the value of N, the number of samples of H(eio
over O::s: co < 2.1t'. Since N is also the period of {hen)}, we must be careful not
to choose N too small. We now illustrate the considerations involved in
choosing the value of N. Let {hen)} be a causal finite-duration sequence
containing M samples. Infinite-duration and noncausal sequences are discussed
later. As shown in Fig. 4.2, if M ::s: N, then

hen) =

{~(n)

forO::s:n::S:M-l
for M::S:n::s:N-l

(4.4)

Two points are noteworthy about this result when M ::s: N. The first is that
the finite-duration sequence h (n), for O::s: n ::s: M - 1, can be recovered
uniquely from the first period of {hen)}. The second point is that the excess
number of points in one period of {h(n)}, or those for kM::s: n ::s:
keN -1) for any integer k, are equal to zero.

{h(n)} (M = 5)

..... rIIr
-10-1-2345

000

{h(n)} (N=5)

-1 0 1 2 3 4 5
{h(n)} (N=8)

.. Irllt ... tlllt tII ~

~3

{ii(n)}

4 567,8

(N = 4)

I I I I

1 1 11I 1I I I I 1 I 1
n

FIGURE 4-2
Relationship of duration of sequence M and number of sample
points in spectrum N. When N <
M, a time-aliasing elIect occurs.
Underbar indicates sequence produced by inverse OFf.

TIffi DISCRETE FOURIER TRANSFORM

121

If M > N, however, the duration of {h(n)} is longer than the period of

{h(n)}. An overlap then occurs when the periodic extension is formed, as

shown in Fig. 4.2. This overlap results in an error when we try to retrieve
{h(n)} from the first period of {h(n)}. In this case, time-aliasing is said to
occur. To prevent time-aliasing, we must sample H(ejOO) over 0:5 OJ < 23f at
least as many times as there are elements in {h(n)}, or N? M times. Clearly,
this poses a difficulty when trying to apply this result to infinite-duration
sequences. These problems will be considered below.

4.3 COMPUTING mE DISCRETE

FOURIER TRANSFORM FROM mE
DISCRETE-TIME SEQUENCE
We now derive the relationship of the discrete Fourier transform to the
discrete-time sequence. First, we note that since {h(n)} is a periodic sequence,
it can be expressed as the following Fourier series expansion
h(n)

=..!.
N

H(k) ej'br1cnIN

(4.5)

k=-oo

where {ii(k)} are a set of complex-valued coefficients to be determined, and

the lIN factor has been accepted as the convention. Because the discrete-time
complex exponential sequence {ei'brknIN} is itself periodic in k, with period N,
there are only a finite number of distinct sequences, equal to ei'brknlN, for
0:5 k s N - 1. To show this, let 0:5 k :5 N - 1 and consider a frequency
outside this range, say ei'br(k+N)nIN. Then
ej'br(k+N)IIIN = ei'br1cnIN ei27rn = ej 'br1cnIN
(4.6)
which is identical to the sequence ei'brknIN. Hence, only a finite set of
frequencies need to be included in the above summation, and {h(n)} can be
expressed in terms of the following N frequency components
1 N-l
h(n) = H(k) ej'brnklN
N k=O

2:

(4.7)

To determine the values of H(k), for 0 s k:5 N -1, we multiply both sides of
Eq. (4.7) by e-i'brnrIN and sum over OSn:5N -1, to get

2: h(n)~-i27rnrIN=_1 2: 2: H(k)ei27rnkIN

N-l

n~

N-l [N-l

Nn~

k~

e-j'brnrIN

(4.8)

Changing the order of summations and noting that R(k) is not a function