Digital-Signal-Processing-DSP-Using-MATLAB-Wavelets.pdf

Loading ...
Global Do...
News & Politics
DIGITAL SIGNAL PROCESSINGDIGITAL SIGNAL 
PROCESSING
Using MATLAB®
and Wavelets
Michael Weeks
Shelving: Engineering / Computer Science 
ISBN: 0-9778582-0-0
Level: Intermediate to Advanced 
U.S. $69.95 / Canada $85.50
All trademarks and service marks are the property of their respective owners.
Cover design: Tyler Creative
INFINITY SCIENCE PRESS
11 Leavitt Street
Hingham, MA 02043
(781) 740-4487
(781) 740-1677 FAX
info@infi nitysciencepress.com
www.infi nitysciencepress.com
Using MATLAB® and WaveletsE L E C T R I C A L  E N G I N E E R I N G  S E R I E S
DIGITAL SIGNAL PROCESSING
Using MATLAB® and Wavelets
Michael Weeks
Although DSP has long been considered an EE topic, recent developments have also generated sig-
nifi cant interest from the computer science community. DSP applications in the consumer market, 
such as bioinformatics, the MP3 audio format, and MPEG-based cable/satellite television have fueled 
a desire to understand this technology outside of hardware circles.
Designed for upper division engineering and computer science students as well as practicing engineers, 
Digital Signal Processing Using MATLAB and Wavelets emphasizes the practical applications of signal 
processing. Over 100 MATLAB examples and wavelet techniques provide the latest applications of DSP, 
including image processing, games, fi lters, transforms, networking, parallel processing, and sound. 
The book also provides the mathematical processes and techniques needed to ensure an understand-
ing of DSP theory. Designed to be incremental in diffi culty, the book will benefi t readers who are 
unfamiliar with complex mathematical topics or those limited in programming experience. Beginning 
with an introduction to MATLAB programming, it moves through fi lters, sinusoids, sampling, the Fou-
rier transform, the z-transform and other key topics. An entire chapter is dedicated to the discussion 
of wavelets and their applications. A CD-ROM (platform independent) accompanies the book and 
contains source code, projects for each chapter, and the fi gures contained in the book.
FEATURES:
■ Contains over 100 short examples in MATLAB used 
throughout the book
■ Includes an entire chapter on the wavelet transform
■ Designed for the reader who does not have extensive 
math and programming experience
■ Accompanied by a CD-ROM containing MATLAB 
examples, source code, projects, and fi gures from 
the book
■ Contains modern applications of DSP and MATLAB 
project ideas
ABOUT THE AUTHOR:
Michael Weeks is an associate professor at Georgia State University where he teaches courses in Digital Signal Processing. 
He holds a PhD in computer engineering from the University of Louisiana at Lafayette and has authored or co-authored 
numerous journal and conference papers.
WEEKS
BRIEF TABLE OF CONTENTS:
1. Introduction. 2. MATLAB. 
3. Filters. 4. Sinusoids. 5. Sampling. 
6. The Fourier Transform. 7. The 
Number e. 8. The z-Transform. 
9. The Wavelet Transform 
10. Applications. Appendix A. 
Constants and Variables. 
B. Equations. C. DSP Project Ideas. 
D. About the CD. Answers. 
Glossary. Index.
weeks_DSP.indd   1
8/11/06   1:15:29 PM
Digital Signal Processing
Using MATLAB r© and Wavelets
License, Disclaimer of Liability, and Limited Warranty
The CD-ROM that accompanies this book may only be used on a single PC. This
license does not permit its use on the Internet or on a network (of any kind). By
purchasing or using this book/CD-ROM package (the “Work”), you agree that this
license grants permission to use the products contained herein, but does not give
you the right of ownership to any of the textual content in the book or ownership to
any of the information or products contained on the CD-ROM. Use of third party
software contained herein is limited to and subject to licensing terms for the respec-
tive products, and permission must be obtained from the publisher or the owner of
the software in order to reproduce or network any portion of the textual material
or software (in any media) that is contained in the Work.
Infinity Science Press LLC (“ISP” or “the Publisher”) and anyone involved
in the creation, writing or production of the accompanying algorithms, code, or
computer programs (“the software”) or any of the third party software contained on
the CD-ROM or any of the textual material in the book, cannot and do not warrant
the performance or results that might be obtained by using the software or contents
of the book. The authors, developers, and the publisher have used their best ef-
forts to insure the accuracy and functionality of the textual material and programs
contained in this package; we, however, make no warranty of any kind, express or
implied, regarding the performance of these contents or programs. The Work is sold
“as is” without warranty (except for defective materials used in manufacturing the
disc or due to faulty workmanship);
The authors, developers, and the publisher of any third party software, and anyone
involved in the composition, production, and manufacturing of this work will not be
liable for damages of any kind arising out of the use of (or the inability to use) the
algorithms, source code, computer programs, or textual material contained in this
publication. This includes, but is not limited to, loss of revenue or profit, or other
incidental, physical, or consequential damages arising out of the use of this Work.
The sole remedy in the event of a claim of any kind is expressly limited to replace-
ment of the book and/or the CD-ROM, and only at the discretion of the Publisher.
The use of “implied warranty” and certain “exclusions” vary from state to state,
and might not apply to the purchaser of this product.
Digital Signal Processing
Using MATLAB r© and Wavelets
Michael Weeks
Georgia State University
Infinity Science Press LLC
Hingham, Massachusetts
Copyright 2007 by Infinity Science Press LLC
All rights reserved.
This publication, portions of it, or any accompanying software may not be reproduced in any way, stored in
a retrieval system of any type, or transmitted by any means or media, electronic or mechanical, including,
but not limited to, photocopy, recording, Internet postings or scanning, without prior permission in writing
from the publisher.
Publisher: David F. Pallai
Infinity Science Press LLC
11 Leavitt Street
Hingham, MA 02043
Tel. 877-266-5796 (toll free)
Fax 781-740-1677
info@infinitysciencepress.com
www.infinitysciencepress.com
This book is printed on acid-free paper.
Michael Weeks. Digital Signal Processing Using MATLAB and Wavelets.
ISBN: 0-9778582-0-0
The publisher recognizes and respects all marks used by companies, manufacturers, and developers as a
means to distinguish their products. All brand names and product names mentioned in this book are trade-
marks or service marks of their respective companies. Any omission or misuse (of any kind) of service marks
or trademarks, etc. is not an attempt to infringe on the property of others.
Library of Congress Cataloging-in-Publication Data
Weeks, Michael.
Digital signal processing using MATLAB and Wavelets / Michael Weeks.
p. cm.
Includes index.
ISBN 0-9778582-0-0 (hardcover with cd-rom : alk. paper)
1.
Signal processing–Digital techniques–Mathematics. 2. MATLAB. 3. Wavelets (Mathematics)
I.
Title.
TK5102.9.W433 2006
621.382’2–dc22
2006021318
6 7 8 9 5 4 3 2 1
Our titles are available for adoption, license or bulk purchase by institutions, corporations, etc. For ad-
ditional information, please contact the Customer Service Dept. at 877-266-5796 (toll free).
Requests for replacement of a defective CD-ROM must be accompanied by the original disc, your mail-
ing address, telephone number, date of purchase and purchase price. Please state the nature of the problem,
and send the information to Infinity Science Press, 11 Leavitt Street, Hingham, MA 02043.
The sole obligation of Infinity Science Press to the purchaser is to replace the disc, based on defec-
tive materials or faulty workmanship, but not based on the operation or functionality of the product.
I dedicate this book to my wife Sophie. Je t’aimerai pour toujours.
Contents
Preface
xxi
1 Introduction
1
1.1 Numbers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
1.1.1 Why Do We Use a Base 10 Number System? . . . . . . . . .
2
1.1.2 Why Do Computers Use Binary?
. . . . . . . . . . . . . . .
2
1.1.3 Why Do Programmers Sometimes Use Base 16
(Hexadecimal)? . . . . . . . . . . . . . . . . . . . . . . . . . .
3
1.1.4 Other Number Concepts . . . . . . . . . . . . . . . . . . . . .
4
1.1.5 Complex Numbers . . . . . . . . . . . . . . . . . . . . . . . .
5
1.2 What Is a Signal?
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
10
1.3 Analog Versus Digital
. . . . . . . . . . . . . . . . . . . . . . . . . .
14
1.4 What Is a System? . . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
1.5 What Is a Transform? . . . . . . . . . . . . . . . . . . . . . . . . . .
20
1.6 Why Do We Study Sinusoids? . . . . . . . . . . . . . . . . . . . . . .
22
1.7 Sinusoids and Frequency Plots
. . . . . . . . . . . . . . . . . . . . .
24
1.8 Summations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
26
1.9 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
1.10 Review Questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
2 MATLAB
29
2.1 Working with Variables
. . . . . . . . . . . . . . . . . . . . . . . . .
30
2.2 Getting Help and Writing Comments
. . . . . . . . . . . . . . . . .
31
2.3 MATLAB Programming Basics . . . . . . . . . . . . . . . . . . . . .
32
2.3.1
Scalars, Vectors, and Matrices
. . . . . . . . . . . . . . . . .
33
2.3.2 Number Ranges . . . . . . . . . . . . . . . . . . . . . . . . . .
35
2.3.3 Output
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
35
2.3.4 Conditional Statements (if) . . . . . . . . . . . . . . . . . . .
36
2.3.5 Loops . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
39
vii
viii
DSP Using MATLAB and Wavelets
2.3.6 Continuing a Line
. . . . . . . . . . . . . . . . . . . . . . . .
39
2.4 Arithmetic Examples . . . . . . . . . . . . . . . . . . . . . . . . . . .
39
2.5 Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
52
2.6 How NOT to Plot a Sinusoid . . . . . . . . . . . . . . . . . . . . . .
53
2.7 Plotting a Sinusoid . . . . . . . . . . . . . . . . . . . . . . . . . . . .
56
2.8 Plotting Sinusoids a Little at a Time
. . . . . . . . . . . . . . . . .
60
2.9 Calculating Error . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
63
2.10 Sometimes 0 Is Not Exactly 0 . . . . . . . . . . . . . . . . . . . . . .
64
2.10.1 Comparing Numbers with a Tolerance . . . . . . . . . . . . .
65
2.10.2 Rounding and Truncating . . . . . . . . . . . . . . . . . . . .
69
2.11 MATLAB Programming Tips . . . . . . . . . . . . . . . . . . . . . .
70
2.12 MATLAB Programming Exercises
. . . . . . . . . . . . . . . . . . .
71
2.13 Other Useful MATLAB Commands . . . . . . . . . . . . . . . . . . .
81
2.14 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
83
2.15 Review Questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
83
3 Filters
85
3.1 Parts of a Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
89
3.2 FIR Filter Structures . . . . . . . . . . . . . . . . . . . . . . . . . . .
91
3.3 Causality, Linearity, and Time-Invariance
. . . . . . . . . . . . . . .
98
3.4 Multiply Accumulate Cells . . . . . . . . . . . . . . . . . . . . . . . . 103
3.5 Frequency Response of Filters . . . . . . . . . . . . . . . . . . . . . . 104
3.6
IIR Filters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
3.7 Trends of a Simple IIR Filter . . . . . . . . . . . . . . . . . . . . . . 113
3.8 Correlation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
3.9 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
3.10 Review Questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
4 Sinusoids
133
4.1 Review of Geometry and Trigonometry . . . . . . . . . . . . . . . . . 133
4.2 The Number π
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
4.3 Unit Circles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
4.4 Principal Value of the Phase Shift
. . . . . . . . . . . . . . . . . . . 138
4.5 Amplitudes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
4.6 Harmonic Signals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
4.7 Representing a Digital Signal as a Sum of
Sinusoids
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
4.8 Spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152
4.9 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156
Contents
ix
4.10 Review Questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156
5 Sampling
159
5.1 Sampling
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160
5.2 Reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162
5.3 Sampling and High-Frequency Noise . . . . . . . . . . . . . . . . . . 162
5.4 Aliasing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164
5.4.1 Aliasing Example . . . . . . . . . . . . . . . . . . . . . . . . . 165
5.4.2 Folding
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168
5.4.3 Locations of Replications After Sampling
. . . . . . . . . . . 171
5.5 Nyquist Rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175
5.6 Bandpass Sampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176
5.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182
5.8 Review Questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 183
6 The Fourier Transform
187
6.1 Fast Fourier Transform Versus the Discrete Fourier Transform . . . . 190
6.2 The Discrete Fourier Transform . . . . . . . . . . . . . . . . . . . . . 191
6.3 Plotting the Spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . 196
6.4 Zero Padding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202
6.5 DFT Shifting Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . 203
6.6 The Inverse Discrete Fourier Transform . . . . . . . . . . . . . . . . 204
6.7 Forward and Inverse DFT . . . . . . . . . . . . . . . . . . . . . . . . 207
6.8 Leakage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 212
6.9 Harmonics and Fourier Transform . . . . . . . . . . . . . . . . . . . 214
6.10 Sampling Frequency and the Spectrum . . . . . . . . . . . . . . . . . 219
6.11 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 221
6.12 Review Questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 221
7 The Number e
225
7.1 Reviewing Complex Numbers . . . . . . . . . . . . . . . . . . . . . . 225
7.2 Some Interesting Properties of j
. . . . . . . . . . . . . . . . . . . . 228
7.2.1 Rotating Counterclockwise
. . . . . . . . . . . . . . . . . . . 228
7.2.2 Rotating Clockwise . . . . . . . . . . . . . . . . . . . . . . . . 229
7.2.3 Removing j from
√−a . . . . . . . . . . . . . . . . . . . . . . 230
7.3 Where Does e Come from?
. . . . . . . . . . . . . . . . . . . . . . . 230
7.4 Euler’s Formula . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233
7.5 Alternate Form of Euler’s Equation . . . . . . . . . . . . . . . . . . . 235
7.6 Euler’s Inverse Formula
. . . . . . . . . . . . . . . . . . . . . . . . . 236
7.7 Manipulating Vectors
. . . . . . . . . . . . . . . . . . . . . . . . . . 238
x
DSP Using MATLAB and Wavelets
7.7.1 Adding Two Vectors . . . . . . . . . . . . . . . . . . . . . . . 238
7.7.2 Adding Vectors in General . . . . . . . . . . . . . . . . . . . . 239
7.7.3 Adding Rotating Phasors . . . . . . . . . . . . . . . . . . . . 240
7.7.4 Adding Sinusoids of the Same Frequency
. . . . . . . . . . . 240
7.7.5 Multiplying Complex Numbers . . . . . . . . . . . . . . . . . 240
7.8 Adding Rotating Phasors: an Example . . . . . . . . . . . . . . . . . 243
7.9 Multiplying Phasors
. . . . . . . . . . . . . . . . . . . . . . . . . . . 249
7.10 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 250
7.11 Review Questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 250
8 The z-Transform
253
8.1 The z-Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 254
8.2 Replacing Two FIR Filters in Series
. . . . . . . . . . . . . . . . . . 255
8.3 Revisiting Sequential Filter Combination with z . . . . . . . . . . . . 257
8.4 Why Is z−1 the Same as a Delay by One Unit? . . . . . . . . . . . . 259
8.5 What Is z?
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 260
8.6 How the z-Transform Reduces to the Fourier
Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 261
8.7 Powers of −z
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 261
8.8 Showing that x[n] ∗ h[n] ↔ X(z)H(z)
. . . . . . . . . . . . . . . . . 262
8.9 Frequency Response of Filters . . . . . . . . . . . . . . . . . . . . . . 263
8.10 Trends of a Simple IIR Filter, Part II
. . . . . . . . . . . . . . . . . 271
8.11 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 271
8.12 Review Questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 272
9 The Wavelet Transform
275
9.1 The Two-Channel Filter Bank . . . . . . . . . . . . . . . . . . . . . . 277
9.2 Quadrature Mirror Filters and Conjugate
Quadrature Filters . . . . . . . . . . . . . . . . . . . . . . . . . . . . 279
9.3 How the Haar Transform Is a 45-Degree Rotation . . . . . . . . . . . 280
9.3.1 How The Haar Transform Affects a Point’s Radius . . . . . . 281
9.3.2 How The Haar Transform Affects a Point’s Angle . . . . . . . 282
9.4 Daubechies Four-Coefficient Wavelet . . . . . . . . . . . . . . . . . . 284
9.5 Down-Sampling and Up-Sampling
. . . . . . . . . . . . . . . . . . . 288
9.5.1 Example Using Down/Up-Samplers
. . . . . . . . . . . . . . 288
9.5.2 Down-Sampling and Up-Sampling with 2 Coefficients . . . . . 290
9.5.3 Down-Sampling and Up-Sampling with Daubechies 4 . . . . . 292
9.6 Breaking a Signal Into Waves . . . . . . . . . . . . . . . . . . . . . . 295
Contents
xi
9.7 Wavelet Filter Design—Filters with Four
Coefficients . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 309
9.8 Orthonormal Bases . . . . . . . . . . . . . . . . . . . . . . . . . . . . 311
9.9 Multiresolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 314
9.10 Biorthogonal Wavelets . . . . . . . . . . . . . . . . . . . . . . . . . . 320
9.11 Wavelet Transform Theory
. . . . . . . . . . . . . . . . . . . . . . . 324
9.12 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 336
9.13 Review Questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 336
10 Applications
339
10.1 Examples Working with Sound . . . . . . . . . . . . . . . . . . . . . 339
10.2 Examples Working with Images . . . . . . . . . . . . . . . . . . . . . 342
10.3 Performing the 2D Discrete Wavelet Transform on an Image . . . . . 344
10.3.1 2D DWT of a Grayscale Image . . . . . . . . . . . . . . . . . 347
10.3.2 2D DWT of a Color Image
. . . . . . . . . . . . . . . . . . . 348
10.4 The Plus/Minus Transform . . . . . . . . . . . . . . . . . . . . . . . 350
10.5 Doing and Undoing the Discrete Wavelet
Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 352
10.6 Wavelet Transform with Matrices . . . . . . . . . . . . . . . . . . . . 357
10.7 Recursively Solving a Su Doku Puzzle
. . . . . . . . . . . . . . . . . 361
10.8 Converting Decimal to Binary . . . . . . . . . . . . . . . . . . . . . . 369
10.9 Frequency Magnitude Response of Sound . . . . . . . . . . . . . . . 373
10.10 Filter Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 376
10.10.1 Windowing Methods
. . . . . . . . . . . . . . . . . . . . . . 376
10.10.2 Designing an FIR Filter . . . . . . . . . . . . . . . . . . . . . 380
10.11 Compression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 385
10.11.1 Experimenting with Compression . . . . . . . . . . . . . . . 387
10.11.2 Compressing an Image Ourselves
. . . . . . . . . . . . . . . 393
10.12 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 396
10.13 Review Questions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 396
A Constants and Variables Used in This Book
399
A.1 Constants . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 399
A.2 Variables
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 399
A.3 Symbols Common in DSP Literature . . . . . . . . . . . . . . . . . . 401
B Equations
403
B.1 Euler’s Formula . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403
B.2 Trigonometric Identities and Other Math Notes . . . . . . . . . . . . 404
B.3 Sampling
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 405
xii
DSP Using MATLAB and Wavelets
B.4 Fourier Transform (FT) . . . . . . . . . . . . . . . . . . . . . . . . . 405
B.5 Convolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 407
B.6 Statistics
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 407
B.7 Wavelet Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . 409
B.8 z-Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 409
C DSP Project Ideas
411
D About the CD-ROM
415
E Answers to Selected Review Questions
417
F Glossary
439
Bibliography
445
Index
449
List of Figures
1.1 An example vector. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
1.2 Calculating θ = arctan(b/a) leads to a problem when a is negative. .
8
1.3 A sound signal with a tape analog. . . . . . . . . . . . . . . . . . . .
15
1.4 Sampling a continuous signal. . . . . . . . . . . . . . . . . . . . . . .
17
1.5 Three ways of viewing a signal. . . . . . . . . . . . . . . . . . . . . .
19
1.6 An example system.
. . . . . . . . . . . . . . . . . . . . . . . . . . .
20
1.7 Three glasses of water. . . . . . . . . . . . . . . . . . . . . . . . . . .
21
2.1 A 200 Hz sinusoid produced by example MATLAB code. . . . . . . .
57
2.2 Using the “plotsinusoids” function. . . . . . . . . . . . . . . . . . . .
59
2.3 This signal repeats itself every second. . . . . . . . . . . . . . . . . .
61
2.4 A close-up view of two sinusoids from 0.9 to 1 second. . . . . . . . .
62
3.1 An example signal, filtered.
. . . . . . . . . . . . . . . . . . . . . . .
86
3.2 The frequency content of the example signal, and low/highpass filters. 87
3.3 A digital signal, delayed, appears as a time-shifted version of itself. .
88
3.4 An adder with two signals as inputs. . . . . . . . . . . . . . . . . . .
89
3.5 A multiplier with two signals as inputs.
. . . . . . . . . . . . . . . .
89
3.6 An example FIR filter with coefficients 0.5 and -0.5.
. . . . . . . . .
90
3.7 Signal y is a delayed version of x. . . . . . . . . . . . . . . . . . . . .
91
3.8 FIR filter with coefficients {0.5, 0.5}. . . . . . . . . . . . . . . . . . .
91
3.9 An example FIR filter with coefficients 0.6 and 0.2. . . . . . . . . . .
92
3.10 General form of the FIR filter.
. . . . . . . . . . . . . . . . . . . . .
94
3.11 An example FIR filter. . . . . . . . . . . . . . . . . . . . . . . . . . .
97
3.12 A representation of an FIR filter. . . . . . . . . . . . . . . . . . . . .
97
3.13 Linear condition 1: scaling property. . . . . . . . . . . . . . . . . . . 100
3.14 Linear condition 2: additive property.
. . . . . . . . . . . . . . . . . 100
3.15 Multiply accumulate cell.
. . . . . . . . . . . . . . . . . . . . . . . . 103
3.16 Multiply accumulate cells as a filter.
. . . . . . . . . . . . . . . . . . 104
3.17 Frequency magnitude response for a lowpass filter.
. . . . . . . . . . 105
xiii
xiv
DSP Using MATLAB and Wavelets
3.18 Frequency magnitude response for a highpass filter. . . . . . . . . . . 105
3.19 Passband, transition band, and stopband, shown with ripples. . . . . 107
3.20 Frequency magnitude response for a bandpass filter.
. . . . . . . . . 108
3.21 Frequency magnitude response for a bandstop filter.
. . . . . . . . . 108
3.22 A notch filter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
3.23 Frequency magnitude response for a bandpass filter with two passbands.109
3.24 A filter with feed-back. . . . . . . . . . . . . . . . . . . . . . . . . . . 111
3.25 Another filter with feed-back. . . . . . . . . . . . . . . . . . . . . . . 112
3.26 A third filter with feed-back.
. . . . . . . . . . . . . . . . . . . . . . 113
3.27 General form of the IIR filter. . . . . . . . . . . . . . . . . . . . . . . 114
3.28 A simple IIR filter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
3.29 Output from a simple IIR filter. . . . . . . . . . . . . . . . . . . . . . 116
3.30 Two rectangles of different size (scale) and rotation.
. . . . . . . . . 126
3.31 Two rectangles represented as the distance from their centers to their
edges.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
3.32 A rectangle and a triangle.
. . . . . . . . . . . . . . . . . . . . . . . 128
3.33 A rectangle and triangle represented as the distance from their centers
to their edges. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
4.1 A right triangle.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
4.2 An example circle.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
4.3 An angle specified in radians. . . . . . . . . . . . . . . . . . . . . . . 136
4.4 An angle specified in radians. . . . . . . . . . . . . . . . . . . . . . . 137
4.5 A 60 Hz sinusoid. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
4.6 A vector of −a at angle θ = a at angle (θ + π). . . . . . . . . . . . . 139
4.7 A harmonic signal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
4.8 A short signal.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144
4.9 The short signal, repeated.
. . . . . . . . . . . . . . . . . . . . . . . 145
4.10 A digital signal (top) and its sum of sinusoids representation (bottom).146
4.11 The first four sinusoids in the composite signal. . . . . . . . . . . . . 147
4.12 The last four sinusoids in the composite signal.
. . . . . . . . . . . . 148
4.13 Frequency magnitude spectrum and phase angles. . . . . . . . . . . . 153
4.14 Spectrum plot: magnitude of x(t) = 2 + 2 cos(2π(200)t). . . . . . . . 154
4.15 Spectrum plot: magnitudes. . . . . . . . . . . . . . . . . . . . . . . . 155
4.16 Spectrum plot: phase angles.
. . . . . . . . . . . . . . . . . . . . . . 155
5.1 Sampling a noisy signal. . . . . . . . . . . . . . . . . . . . . . . . . . 163
5.2 Aliasing demonstrated. . . . . . . . . . . . . . . . . . . . . . . . . . . 167
5.3
cos(−2π10t − π/3) and cos(2π10t + π/3) produce the same result.
. 170
5.4 Replications for f1 = 1 Hz and fs = 4 samples/second. . . . . . . . . 172
List of Figures
xv
5.5 Replications for f1 = 2.5 Hz and fs = 8 samples/second. . . . . . . . 173
5.6 Replications for f1 = 3 Hz and fs = 4 samples/second. . . . . . . . . 173
5.7 Replications for f1 = 5 Hz and fs = 4 samples/second. . . . . . . . . 174
5.8 Replications for f1 = 3 or 5 Hz and fs = 4 samples/second.
. . . . . 175
5.9 Example signal sampled at 2081 Hz. . . . . . . . . . . . . . . . . . . 177
5.10 Example signal sampled at 100 Hz. . . . . . . . . . . . . . . . . . . . 177
5.11 Example signal sampled at 103 Hz. . . . . . . . . . . . . . . . . . . . 179
5.12 Example signal sampled at 105 Hz. . . . . . . . . . . . . . . . . . . . 181
5.13 Example signal sampled at 110 Hz. . . . . . . . . . . . . . . . . . . . 181
6.1 A person vocalizing the “ee” sound.
. . . . . . . . . . . . . . . . . . 188
6.2 J.S. Bach’s Adagio from Toccata and Fuge in C—frequency magni-
tude response. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 189
6.3 A sustained note from a flute. . . . . . . . . . . . . . . . . . . . . . . 190
6.4 Comparing Nlog2(N) (line) versus N
2 (asterisks).
. . . . . . . . . . 191
6.5 Spectrum for an example signal.
. . . . . . . . . . . . . . . . . . . . 198
6.6
Improved spectrum for an example signal. . . . . . . . . . . . . . . . 200
6.7 Example output of DFT-shift program.
. . . . . . . . . . . . . . . . 205
6.8 Frequency content appears at exact analysis frequencies. . . . . . . . 213
6.9 Frequency content appears spread out over analysis frequencies. . . . 213
6.10 Approximating a triangle wave with sinusoids. . . . . . . . . . . . . . 215
6.11 Approximating a saw-tooth wave with sinusoids.
. . . . . . . . . . . 217
6.12 Approximating a square wave with sinusoids.
. . . . . . . . . . . . . 218
6.13 Approximating a saw-tooth square wave with sinusoids.
. . . . . . . 218
6.14 Frequency response of a lowpass filter. . . . . . . . . . . . . . . . . . 220
6.15 Frequency response of a highpass filter.
. . . . . . . . . . . . . . . . 220
7.1 A complex number can be shown as a point or a 2D vector. . . . . . 226
7.2 A vector forms a right triangle with the x-axis.
. . . . . . . . . . . . 226
7.3 A rotating vector.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . 227
7.4 A vector rotates clockwise by −π/2 when multiplied by −j. . . . . . 230
7.5 Simple versus compounded interest.
. . . . . . . . . . . . . . . . . . 232
7.6 Phasors and their complex conjugates. . . . . . . . . . . . . . . . . . 236
7.7 Graph of x(0) where x(t) = 3ejπ/6ej2π1000t.
. . . . . . . . . . . . . . 237
7.8 Two example vectors.
. . . . . . . . . . . . . . . . . . . . . . . . . . 241
7.9 Adding two example vectors.
. . . . . . . . . . . . . . . . . . . . . . 241
7.10 Adding and multiplying two example vectors. . . . . . . . . . . . . . 242
7.11 Two sinusoids of the same frequency added point-for-point and ana-
lytically. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 245
7.12 A graphic representation of adding 2 phasors of the same frequency.
247
xvi
DSP Using MATLAB and Wavelets
7.13 A graphic representation of adding 2 sinusoids of the same frequency. 248
8.1 FIR filters in series can be combined. . . . . . . . . . . . . . . . . . . 256
8.2 Two trivial FIR filters. . . . . . . . . . . . . . . . . . . . . . . . . . . 259
8.3 Two trivial FIR filters, reduced. . . . . . . . . . . . . . . . . . . . . . 260
8.4 Example plot of zeros and poles.
. . . . . . . . . . . . . . . . . . . . 270
9.1 Analysis filters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 276
9.2 Synthesis filters.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 276
9.3 A two-channel filter bank. . . . . . . . . . . . . . . . . . . . . . . . . 277
9.4 A quadrature mirror filter for the Haar transform.
. . . . . . . . . . 280
9.5 A conjugate quadrature filter for the Haar transform. . . . . . . . . . 280
9.6 A two-channel filter bank with four coefficients. . . . . . . . . . . . . 284
9.7 Different ways to indicate down-samplers and up-samplers.
. . . . . 288
9.8 A simple filter bank demonstrating down/up-sampling. . . . . . . . . 288
9.9 A simple filter bank demonstrating down/up-sampling, reduced.
. . 289
9.10 Tracing input to output of a simple filter bank. . . . . . . . . . . . . 289
9.11 A two-channel filter bank with down/up-samplers.
. . . . . . . . . . 290
9.12 A filter bank with four taps per filter.
. . . . . . . . . . . . . . . . . 292
9.13 Wavelet analysis and reconstruction. . . . . . . . . . . . . . . . . . . 296
9.14 Alternate wavelet reconstruction. . . . . . . . . . . . . . . . . . . . . 296
9.15 Impulse function analyzed with Haar.
. . . . . . . . . . . . . . . . . 298
9.16 Impulse function analyzed with Daubechies-2. . . . . . . . . . . . . . 299
9.17 Impulse function (original) and its reconstruction.
. . . . . . . . . . 300
9.18 Example function broken down into 3 details and approximation. . . 301
9.19 Impulse function in Fourier-domain.
. . . . . . . . . . . . . . . . . . 304
9.20 Impulse function in Wavelet-domain. . . . . . . . . . . . . . . . . . . 305
9.21 Designing a filter bank with four taps each.
. . . . . . . . . . . . . . 310
9.22 Two levels of resolution. . . . . . . . . . . . . . . . . . . . . . . . . . 315
9.23 Biorthogonal wavelet transform. . . . . . . . . . . . . . . . . . . . . . 320
9.24 Biorthogonal wavelet transform. . . . . . . . . . . . . . . . . . . . . . 322
9.25 One channel of the analysis side of a filter bank.
. . . . . . . . . . . 329
9.26 The first graph of the “show wavelet” program. . . . . . . . . . . . . 334
9.27 The second graph of the “show wavelet” program.
. . . . . . . . . . 335
9.28 Analysis filters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 337
9.29 Synthesis filters.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 337
10.1 Two-dimensional DWT on an image. . . . . . . . . . . . . . . . . . . 345
10.2 Results of “DWT undo” for a 1D signal: original and approximation. 354
10.3 Results of “DWT undo” for a 1D signal: details for octaves 1–3.
. . 355
List of Figures
xvii
10.4 Three octaves of the one-dimensional discrete wavelet transform. . . 362
10.5 Example of the “show sound” program.
. . . . . . . . . . . . . . . . 377
10.6 Two similar filters, with and without a gradual transition. . . . . . . 379
10.7 Using nonwindowed filter coefficients. . . . . . . . . . . . . . . . . . . 382
10.8 Using windowed filter coefficients. . . . . . . . . . . . . . . . . . . . . 383
10.9 Windowed filter coefficients and those generated by fir1.
. . . . . . 384
10.10 Alternating flip for the windowed filter coefficients. . . . . . . . . . 386
List of Tables
1.1 Example signals.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
12
5.1 Frequencies detected with bandpass sampling. . . . . . . . . . . . . . 182
6.1 Example DFT calculations.
. . . . . . . . . . . . . . . . . . . . . . . 194
6.2 Example DFT calculations (rectangular form).
. . . . . . . . . . . . 195
6.3 Sinusoids that simplify things for us. . . . . . . . . . . . . . . . . . . 211
7.1 Compound interest on $1000 approximates $1000 ×e.
. . . . . . . . 232
8.1 Convolution of x and h.
. . . . . . . . . . . . . . . . . . . . . . . . . 263
9.1
Impulse function analyzed with Haar.
. . . . . . . . . . . . . . . . . 298
10.1 An example Su Doku puzzle.
. . . . . . . . . . . . . . . . . . . . . . 363
10.2 The example Su Doku puzzle’s solution. . . . . . . . . . . . . . . . . 364
10.3 Converting from decimal to binary. . . . . . . . . . . . . . . . . . . . 370
10.4 Converting from a decimal fraction to fixed-point binary. . . . . . . . 371
10.5 Hexadecimal to binary chart.
. . . . . . . . . . . . . . . . . . . . . . 374
A.1 Greek alphabet. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 402
xix
Preface
Digital signal processing is an important and growing field. There have been
many books written in this area, however, I was motivated to write this manuscript
because none of the existing books address the needs of the computer scientist. This
work attempts to fill this void, and to bridge the disciplines from which this field
originates: mathematics, electrical engineering, physics, and engineering mechanics.
Herein, it is hoped that the reader will find sensible explanations of DSP concepts,
along with the analytical tools on which DSP is based.
DSP: Now and Then
Recently, I was asked if a student who took a course in DSP at another university
should be allowed to transfer the credit to our university. The question was an
interesting one, especially when I noticed that the student had taken the course 10
years earlier. For computer science, 10 years is a long time. But how much has DSP
changed in the last decade? It occurred to me that there were several important
changes, but theoretically, the core information has not changed. The curriculum
of the DSP course in question contained all of the same topics that I covered in my
DSP class that summer, with only two exceptions: MATLAB r© and Wavelets.
MATLAB
MATLAB is one of several tools for working with mathematics. As an experienced
programmer, I was skeptical at first. Why should I learn another programming
language when I could do the work in C/C++? The answer was simple: working
with MATLAB is easier! Yes, there are some instances when one would want to
use another language. A compiled program, such as one written in C++, will run
faster than an interpreted MATLAB program (where each line is translated by the
computer when the program is run). For hardware design, one might want to write
a program in Verilog or VHDL, so that the program can be converted to a circuit. If
you already know a programming language such as C, C++, Java, FORTRAN, etc.,
xxi
xxii
DSP Using MATLAB and Wavelets
you should be able to pick up a MATLAB program and understand what it does. If
you are new to programming, you will find MATLAB to be a forgiving environment
where you can test out commands and correct them as needed.
Wavelets
The wavelet transform is an analysis tool that has a relatively short history. It
was not until 1989 that Stéphane Mallat (with some help from Meyers) published
a revolutionary paper on wavelet theory, tying together related techniques that
existed in several disciplines [1]. Other people have made significant contributions
to this theory, such as Ingrid Daubechies [2], Strang and Nguyen [3], and Coifman
and Wickerhauser [4]. It is a daunting task to write about the people who have
contributed to wavelets, since anyone who is left out could be reading this now!
The wavelet transform is an important tool with many applications, such as
compression. I have no doubt that future generations of DSP teachers will rank it
second only to the Fourier transform in terms of importance.
Other Developments
Other recent changes in DSP include embedded systems, the audio format MP3, and
public awareness. Advertising by cellular phone marketers, which tries to explain
to a nontechnical audience that digital signal processing is better than analog, is an
example of the growing public awareness. Actually, the ads do not try to say how
digital is better than analog, but they do point out problems with wireless analog
phones.
Acknowledgments
This book was only possible with the support of many people.
I would like to
thank Dr. Rammohan Ragade from the University of Louisville, who guided my
research when I first started out. I would also like to thank Dr. Magdy Bayoumi
at the University of Louisiana, who taught me much about research as well as how
to present research results to an audience. Dr. Marty Fraser, who recently retired
from Georgia State University, was a great mentor in my first years in academia.
My students over the years have been quite helpful in pointing out any problems
within the text. I would especially like to thank Evelyn Brannock and Ferrol Black-
mon for their help reviewing the material. I would also like to acknowledge Drs.
Preface
xxiii
Kim King and Raj Sunderraman for their help. Finally, I could not have writ-
ten this book without the support and understanding of my wife.
-M.C.W., Atlanta, Georgia, 2006
Chapter 1
Introduction
Digital Signal Processing (DSP) is an important field of study that has come
about due to advances in communication theory, digital (computer) technology, and
consumer devices. There is always a driving need to make things better and DSP
provides many techniques for doing this. For example, people enjoy music and like
to download new songs. However, with slow Internet connection speeds (typically
56 kilobits per second for a dial-up modem), downloading a song could take hours.
With MP3 compression software, though, the size of the song is reduced by as much
as 90%, and can be downloaded in a matter of minutes. The MP3 version of the
song is not the same as the original, but is a “good enough” approximation that
most users cannot distinguish from the original. How is this possible? First, it is
important to know about the original song (a signal), and how it is represented
digitally. This knowledge leads to an algorithm to remove data that the user will
not miss. All of this is part of Digital Signal Processing.
To understand how to manipulate (process) a signal, we must first know a bit
about the values that make up a signal.
1.1 Numbers
Consider our concepts of numbers.
In the ancient world, people used numbers
to count things. In fact, the letter “A” was originally written upside down, and
represented an ox [5]. (It looks a bit like an ox’s head if you write it upside down.)
Ratios soon followed, since some things were more valuable than others. If you were
an ancient pig farmer, you might want a loaf of bread, but would you be willing to
trade a pig for a loaf of bread? Maybe you would be willing to trade a pig for a
loaf of bread and three ducks. The point is, ratios developed as a way to compare
dissimilar things. With ratios, come fractions. A duck may be worth three loaves of
1
2
DSP Using MATLAB and Wavelets
bread, and if you traded two ducks you would expect six loaves of bread in return.
This could work the other way, too. If you had a roasted duck, you might divide it
into 3 equal sections, and trade one of the sections for a bread loaf. Ratios lead the
way to fractions, and our use of the decimal point to separate the whole part of the
number from the fractional part.
Zero is a useful number whose invention is often credited to the Arabs, though
there is much debate on this point. Its origins may not ever be conclusive. It is one
of the symbols we use in our base 10 system, along with the numerals 1 through 9.
(Imagine if we used Roman numerals instead!) The concept of zero must have been
revolutionary in its time, because it is an abstract idea. A person can see 1 duck or
3 loaves of bread, but how do you see 0 of something? Still, it is useful, even if just
a counting system is used. We also use it as a placeholder, to differentiate 10 from
1 or 100.
Ten is a significant number that we use as the basis of our decimal number
system. There is no symbol for ten like there is for 0–9 (at least, not in the decimal
number system). Instead, we use a combination of symbols, i.e., 10. We understand
that the 1 is for the ten’s digit, and the 0 is the one’s digit, so that placement is
important. In any number system, we have this idea of placement, where the digits
on the left are greater in magnitude than the digits on the right. This is also true
in number systems that computers use, i.e., binary.
1.1.1 Why Do We Use a Base 10 Number System?
It’s most likely that we use a base 10 number system because humans have 10 fingers
and used them to count with, as children still do today.
1.1.2 Why Do Computers Use Binary?
Computers use base 2, binary, internally. This number system works well with the
computer’s internal parts; namely, transistors. A transistor works as a switch, where
electricity flows easily when a certain value is applied to the gate, and the flow of
electricity is greatly impeded when another value is applied. By viewing the two
values that work with these transistors as a logical 0 (ground) and a logical 1 (e.g.,
3.3 volts, depending on the system), we have a number system of two values. Either
a value is “on,” or it is “off.” One of the strengths of using a binary system is that
correction can take place. For example, if a binary value appears as 0.3 volts, a
person would conclude that this value is supposed to be 0.0 volts, or logic 0. To the
computer, such a value given to a digital logical circuit would be close enough to 0.0
volts to work the same, but would be passed on by the circuit as a corrected value.
Introduction
3
1.1.3 Why Do Programmers Sometimes Use Base 16
(Hexadecimal)?
Computers use binary, which is great for a computer, but difficult for a human. For
example, which is easier for you to remember, 0011001000010101 or 3215? Actually,
these two numbers are the same, the first is in binary and the second is in hexadec-
imal. As you can see, hexadecimal provides an easier way for people to work with
computers, and it translates to and from binary very easily. In fact, one hexadecimal
digit represents four bits. With a group of four bits, the only possibilities are: 0000,
0001, 0010, 0011, 0100, 0101, 0110, 0111, 1000, 1001, 1010, 1011, 1100, 1101, 1110,
and 1111. These 16 possible combinations of 4 bits map to the numbers 0, 1, 2, 3, 4,
5, 6, 7, 8, 9. At this point, we have a problem, since we would be tempted to make
the next number “10,” as we would in decimal. The problem with this is that it is
ambiguous: would “210” mean 2, 1, 0, or would it mean 2, 10? Since we have 16
different values in a group of four binary digits, we need 16 different symbols, one
symbol for each value. This is why characters from the alphabet are used next: A,
B, C, D, E, and F.
It is inefficient to use 1-bit registers, so computers group bits together into a
word. The word size is dependent on the architecture (e.g., the “64” in the Nin-
tendo 64 stands for the 64-bit word size of its processor). Typically, the word size
is a multiple of a byte (8 bits), and hexadecimal numbers work nicely with such
machines. For example, a machine with a 1-byte word size would work with data
sizes of 2 hexadecimal digits at a time. A 16-bit word size means 4 hexadecimal
digits are used.
Here is an example demonstrating word size. Suppose we have two multiplica-
tions,
9 × 3
4632 × 9187.
We can do both of the above multiplications, but most people have an immediate
answer for the one on the left, but need a minute or two for the one on the right.
Why? While they might have the single-digit multiplication answers memorized,
they use an algorithm for the multiple digit multiplications, i.e., multiply the right-
most digits, write the result, write the carry, move one digit to the left, for the
“top” number, and repeat the process. In a sense, we have a 1-digit word size when
performing this calculation. Similarly, a computer will multiply (or add, subtract,
etc.) in terms of its word size. For example, suppose that a computer has a word
size of 8 bits. If it needed to increment a 16-bit number, it would add one to the
low 8 bits, then add the carry to the high 8 bits.
4
DSP Using MATLAB and Wavelets
1.1.4 Other Number Concepts
Negative numbers are another useful concept that comes about due to borrowing. It
is not difficult to imagine someone owing you a quantity of some item. How would
you represent this? A negative amount is what we use today, e.g., in a checking
account, to denote that not only do you have zero of something, but that you owe
a quantity to someone else.
Infinity is another strange but useful idea. With our number system, no matter
what extremely large number someone comes up with, you could still add 1 to it.
So we use infinity to fill in for “the largest possible number.” We can talk about all
numbers by including negative infinity and infinity as the limits of our range. This
is often seen in formulas as a way of covering all numbers.
Most of the discussion will be centered on fixed-point numbers. This means
exactly what the name implies—the radix point is fixed. A radix point is the name
given to the symbol separating the whole part from the fractional part—a period in
the U.S., while Europeans use a comma. It is not exactly a “decimal point,” since
this term implies that base 10 (decimal) is being used. A fixed-point number can be
converted to decimal (and back) in the same way that a binary number is converted.
Recall that to convert a binary number to decimal, one would multiply each bit
in sequence by a multiple of 2, working from the right to the left. Thus, for the
binary number 01101001, the decimal equivalent is 0 × 27 + 1 × 26 + 1 × 25 + 0 ×
24 + 1 × 23 + 0 × 22 + 0 × 21 + 1 × 20 = 1 × 64 + 1 × 32 + 1 × 8 + 1 × 1 = 105.
For a fixed-point number, this idea holds, except that the bit to the left of the
radix point corresponds to 20. One can start here and work left, then start again
at the radix point and work right. For example, converting the fixed-point number
01010.111 would be 0×24+1×23+0×22+1×21+0×20+1×2−1+1×2−2+1×2−3 =
1 × 8 + 1 × 2 + 1 × (1/2) + 1 × (1/4) + 1 × (1/8) = 10.875 in decimal.
To find a decimal number in binary, we separate the whole part from the frac-
tional part. We take the whole part, divide it by 2, and keep track of the result and
its remainder. Then we repeat this with the result until the result becomes zero.
Read the remainders back, from bottom to top, and that is our binary number. For
example, suppose we want to convert the decimal number 4 to binary:
4/2 = 2 remainder 0
2/2 = 1 remainder 0
1/2 = 0 remainder 1.
Therefore, decimal 4 equals 100 in binary. We can precede our answer by 0 to
avoid confusion with a negative binary number, i.e., 0100 binary.
Introduction
5
For a fractional decimal number, multiply by 2 and keep the whole part. Then
repeat with the remaining fractional part until it becomes zero (though it is possible
for it to repeat forever, just as 1/3 does in decimal). When finished, read the whole
parts back from the top down. For example, say we want to convert .375 to binary:
.375 × 2 = 0 plus .75
.75 × 2 = 1 plus .5
.5 × 2 = 1 plus 0.
Our answer for this is .011 binary. We can put the answers together and conclude
that 4.375 equals 0100.011 in binary.
1.1.5 Complex Numbers
Complex numbers are important to digital signal processing. For example, functions
like the Fourier transform (fft and ifft) return their results as complex numbers.
This topic is covered in detail later in this book, but for the moment you can think
of the Fourier transform as a function that converts data to an interesting form.
Complex numbers provide a convenient way to store two pieces of information,
either x and y coordinates, or a magnitude (length of a 2D vector) and an angle
(how much to rotate the vector). This information has a physical interpretation in
some contexts, such as corresponding to the magnitude and phase angle of a sinusoid
for a given frequency, when returned from the Fourier transform.
Most likely, you first ran into complex numbers in a math class, factoring roots
with the quadratic formula. That is, an equation such as
2x2 + 4x− 30 = 0
could be factored as (x− root1)(x− root2), where
root1,2 =
−b±
√
b2 − 4ac
2a
.
If the numbers work out nicely, then the quadratic formula results in a couple
of real roots. In this case, it would be 3,−5, to make (x− 3)(x+5) = 0. A problem
arises when b2 − 4ac is negative because the square root of a negative quantity does
not exist. The square root operator is supposed to return the positive root. How
can we take a positive number, multiply it by itself, and end up with a negative
number? For every “real” value, even the negative ones, the square is positive. To
deal with this, we have the imaginary quantity j. (Some people prefer i instead.)
6
DSP Using MATLAB and Wavelets
This imaginary quantity is defined:
j =
√
−1.
Numbers that have j are called “imaginary” numbers, as well as complex num-
bers. A complex number has a real part and an imaginary part, in the form
real + imaginary × j. For example, x = 3 + 4j. To specify the real or imaginary
part, we use functions of the same name, i.e., Real(x) = 3, and Imaginary(x) = 4,
or simply <(x) = 3, and =(x) = 4.
The best way to think of complex numbers is as a two-dimensional generalization
of numbers; to think of real numbers as having an imaginary component of zero.
Thus, real numbers lie on the x-axis.
Complex numbers can be used to hold two pieces of information. For example,
you probably remember converting numbers with polar coordinates to Cartesian1
ones, which are two ways of expressing the same thing. A complex number can
represent a number in polar coordinates. The best way to show this is graphically.
When we want to draw a point in two-dimensional space, we need to specify an X
coordinate and a Y coordinate. This set of two numbers tells us where the point is
located. Polar coordinates do this too, but the point is specified by a length and an
angle. That is, the length specifies a line segment starting at point 0 and going to
the right, and the angle says how much to rotate the line segment counterclockwise.
In a similar way, a complex number represents a point on this plane, where the
real part and imaginary part combine to specify the point. Figure 1.1 shows a 2D
“vector” r units long, rotated at an angle of θ.
θ
r
Figure 1.1: An example vector.
1Named for Rene Descartes, the same man famous for the saying “cogito ergo sum,” or “I think
therefore I am.”
Introduction
7
To convert from polar coordinates to complex Cartesian ones [6]:
x = r cos(θ)
y = r sin(θ).
The MATLAB r© function below converts polar coordinates to complex coordi-
nates. MATLAB will be explained in Chapter 2, “MATLAB,” but anyone with
some programming experience should be able to get the idea of the code below.
% Convert from magnitude,angle form to x+jy form
%
% usage:
[X] = polar2complex(magnitudes, angles)
function [X] = polar2complex(mag, angles)
% Find x and y coordinates as a complex number
X = mag.*cos(angles) + j*mag.*sin(angles);
To convert from Cartesian coordinates to polar ones:
r = |x+ jy| =
√
x2 + y2
θ = arctan(y/x).
There is a problem with the equation used to convert to polar coordinates.
If both real and imaginary parts (x and y) are negative, then the signs cancel,
and the arctan function returns the same value as if they both were positive. In
other words, converting 2 + j2 to polar coordinates gives the exact same answer as
converting −2− j2. In a similar fashion, the arctan function returns the same value
for arctan(y/−x) as it does for arctan(−y/x). To fix this problem, examine the real
part (x), and if it is negative, add π (or 180◦) to the angle. Note that subtracting
π also works, since the angle is off by π, and the functions are 2π periodic. That is,
the cosine and sine functions return the same results for an argument of θ+π− 2π,
as they do for an argument of θ + π, and θ + π − 2π = θ − π.
Figure 1.2 demonstrates this. A vector is shown in all four quadrants. The
radius r equals
√
x2 + y2, and since x is always positive or negative a, and y is
always positive or negative b, all four values of r are the same, where a and b are
8
DSP Using MATLAB and Wavelets
the distances from the origin along the real and imaginary axes, respectively. Let
a = 3 and b = 4:
.
□□
□□




real axis
imaginary axis
j
θ
z  = a + jb
r1
1






real axis
imaginary axis
real axis
imaginary axis
j
r
j
z  = a − jb
z  = −a − jb
r3
4
3
4
3
real axis
imaginary axis
j
z  = −a + jb
r
2
2
2
θ
θ
θ 4
1
Figure 1.2: Calculating θ = arctan(b/a) leads to a problem when a is negative.
r1 =
√
a2 + b2 =
√
32 + 42 = 5
r2 =
√
(−a)2 + (b)2 =
√
(−3)2 + 42 = 5
r3 =
√
(−a)2 + (−b)2 =
√
(−3)2 + (−4)2 = 5
r4 =
√
(a)2 + (−b)2 =
√
(−3)2 + (−4)2 = 5.
Here are the calculated angles:
Introduction
9
θ1 = arctan(b/a) = arctan(4/3) = 0.9273 rad
θ2 = arctan(b/− a) = arctan(−4/3) = −0.9273 rad
θ3 = arctan(−b/− a) = arctan(4/3) = 0.9273 rad
θ4 = arctan(−b/a) = arctan(−4/3) = −0.9273 rad.
Clearly, θ2 and θ3 are not correct since they lie in the second and third quadrants,
respectively. Therefore, their angles should measure between π/2 (≈ 1.57) and π
(≈ 3.14) for θ2 and between −π/2 and −π for θ3. Adding π to θ2 and −π to θ3 fixes
the problem. θ4 is fine, even though the arctan function returns a negative value for
it. Here are the corrected angles:
θ1 = arctan(b/a) = arctan(4/3) = 0.9273 rad
θ2 = arctan(b/− a) + π = −0.9273 + π = 2.2143 rad
θ3 = arctan(−b/− a) − π = 0.9273 − π = −2.2143 rad
θ4 = arctan(−b/a) = arctan(−4/3) = −0.9273 rad.
The function below converts from complex numbers like x+ jy to polar coordi-
nates. This is not as efficient as using abs and angle, but it demonstrates how to
implement the equations.
%
% Convert from complex form (x+jy)
% to polar form (magnitude,angle)
%
% usage:
[r, theta] = complex2polar(X)
%
function [mag, phase] = complex2polar(X)
% Find magnitudes
mag = sqrt(real(X).*real(X) + imag(X).*imag(X));
% Find phase angles
% Note that parameters for tan and atan are in radians
10
DSP Using MATLAB and Wavelets
for i=1:length(X)
if (real(X(i)) > 0)
phase(i) = atan(imag(X(i)) / real(X(i)));
elseif (real(X(i)) < 0)
% Add to +/- pi, depending on quadrant of the point
if (imag(X(i)) < 0) % then we are in quadrant 3
phase(i) = -pi + atan(imag(X(i)) / real(X(i)));
else
% we are in quadrant 2
phase(i) = pi + atan(imag(X(i)) / real(X(i)));
end
else
% If real part is 0, then it lies on the + or - y-axis
% depending on the sign of the imaginary part
phase(i) = sign(imag(X(i)))*pi/2;
end
end
1.2 What Is a Signal?
A signal is a varying phenomenon that can be measured. It is often a physical
quantity that varies with time, though it could vary with another parameter, such
as space. Examples include sound (or more precisely, acoustical pressure), a voltage
(such as the voltage differences produced by a microphone), radar, and images trans-
mitted by a video camera. Temperature is another example of a signal. Measured
every hour, the temperature will fluctuate, typically going from a cold value (in the
early morning) to a warmer one (late morning), to an even warmer one (afternoon),
to a cooler one (evening) and finally a cold value again at night. Often, we must
examine the signal over a period of time. For example, if you are planning to travel
to a distant city, the city’s average temperature may give you a rough idea of what
clothes to pack. But if you look at how the temperature changes over a day, it will
let you know whether or not you need to bring a jacket.
Signals may include error due to limitations of the measuring device, or due to
the environment. For example, a temperature sensor may be affected by wind chill.
At best, signals represented by a computer are good approximations of the original
physical processes.
Some real signals can be measured continuously, such as the temperature. No
matter what time you look at a thermometer, it will give a reading, even if the
time between the readings is arbitrarily small. We can record the temperature at
intervals of every second, every minute, every hour, etc. Once we have recorded these
Introduction
11
measurements, we understand intuitively that the temperature has values between
readings, and we do not know what values these would be. If a cold wind blows, the
temperature goes down, and if the sun shines through the clouds, then it goes up.
For example, suppose we measure the temperature every hour. By doing this, we
are choosing to ignore the temperature for all time except for the hourly readings.
This is an important idea: the signal may vary over time, but when we take periodic
readings of the signal, we are left with only a representation of the signal.
A signal can be thought of as a (continuous or discrete) sequence of (continuous
or discrete) values. That is, a continuous signal may have values at any arbitrary
index value (you can measure the temperature at noon, or, if you like, you can
measure it at 0.0000000003 seconds after noon). A discrete signal, however, has
restrictions on the index, typically that it must be an integer. For example, the
mass of each planet in our solar system could be recorded, numbering the planets
according to their relative positions from the sun. For simplicity, a discrete signal is
assumed to have an integer index, and the relationship between the index and time
(or whatever parameter) must be given. Likewise, the values for the signal can be
with an arbitrary precision (continuous), or with a restricted precision (discrete).
That is, you could record the temperature out to millionths of a degree, or you
could restrict the values to something reasonable like one digit past the decimal.
Discrete does not mean integer, but rather that the values could be stored as a
rational number (an integer divided by another integer). For example, 72.3 degrees
Fahrenheit could be thought of as 723/10. What this implies is that irrational
numbers cannot be stored in a computer, but only approximated. π is a good
example. You might write 3.14 for π, but this is merely an approximation. If you
wrote 3.141592654 to represent π, this is still only an approximation. In fact, you
could write π out to 50 million digits, but it would still be only an approximation!
It is possible to consider a signal whose index is continuous and whose values
are discrete, such as the number of people who are in a building at any given time.
The index (time) may be measured in fractions of a second, while the number of
people is always a whole number. It is also possible to have a signal where the index
is discrete, and the values are continuous; for example, the time of birth of every
person in a city. Person #4 might have been born only 1 microsecond before person
#5, but they technically were not born at the same time. That does not mean that
two people cannot have the exact same birth time, but that we can be as precise as
we want with this time. Table 1.1 gives a few example signals, with continuous as
well as discrete indices and quantities measured.
For the most part, we will concentrate on continuous signals (which have a
continuous index and a continuous value), and discrete signals (with an integer index
and a discrete value). Most signals in nature are continuous, but signals represented
12
DSP Using MATLAB and Wavelets
Table 1.1: Example signals.
index
quantity measured
discrete
continuous
discrete mass (planet number) people in a building (time)
continuous
birth time (person)
temperature (time)
inside a computer are discrete. A discrete signal is often an approximation of a
continuous signal. One notational convention that we adopt here is to use x[n] for a
discrete signal, and x(t) for a continuous signal. This is useful since you are probably
already familiar with seeing the parentheses for mathematical functions, and using
square brackets for arrays in many computer languages.
Therefore, there are 4 kinds of signals:
• A signal can have a continuous value for a continuous index. The “real world”
is full of such signals, but we must approximate them if we want a digi-
tal computer to work with them. Mathematical functions are also continu-
ous/continuous.
• A signal can have a continuous value for a discrete index.
• A signal can have a discrete value for a continuous index.
• A signal can have a discrete value for a discrete index. This is normally the
case in a computer, since it can only deal with numbers that are limited in
range. A computer could calculate the value of a function (e.g., sin(x)), or it
could store a signal in an array (indexed by a positive integer). Technically,
both of these are discrete/discrete signals.
We will
concentrate on the continuous/continuous
signals,
and the
discrete/discrete signals, since these are what we have in the “real world” and the
“computer world,” respectively. We will hereafter refer to these signals as “analog”
and “digital,” respectively.
In the digital realm, a signal is nothing more than an array of numbers.
It
can be thought of in terms of a vector, a one-dimensional matrix. Of course, there
are multidimensional signals, such as images, which are simply two-dimensional
matrices. Thinking of signals as matrices is an important analytical step, since it
allows us to use linear algebra with our signals. The meaning of these numbers
depends on the application as well as the time difference between values. That
Introduction
13
is, one array of numbers could be the changes of acoustic pressure measured at 1
millisecond intervals, or it could be the temperature in degrees centigrade measured
every hour. This essential information is not stored within the signal, but must be
known in order to make sense of the signal. Signals are often studied in terms of time
and amplitude. Amplitude is used as a general way of labeling the units of a signal,
without being limited by the specific signal. When speaking of the amplitude of a
value in a signal, it does not matter if the value is measured in degrees centigrade,
pressure, or voltage.
Usually, in this text, parentheses will be used for analog signals, and square
brackets will be used for digital signals. Also, the variable t will normally be a
continuous variable denoting time, while n will be discrete and will represent sample
number. It can be thought of as an index, or array offset. Therefore, x(t) will
represent a continuous signal, but x[n] will represent a discrete signal. We could let
t have values such as 1.0, 2.345678, and 6.00004, but n would be limited to integer
values like 1, 2, and 6. Theoretically, n could be negative, though we will limit it to
positive integers {0, 1, 2, 3, etc. } unless there is a compelling reason to use negatives,
too. The programming examples present a case where these conventions will not
be followed. In MATLAB, as in other programming languages, the square brackets
and the parentheses have different meanings, and cannot be used interchangeably.
Any reader who understands programming should have no trouble with signals.
Typically, when we talk about a signal, we mean something concrete, such as mySig-
nal in the following C++ program.
#include <iostream>
using namespace std;
int main()
float mySignal[] = 4.0, 4.6, 5.1, 0.6, 6.0 ;
int n, N=5;
// Print the signal values
for (n=0; n<N; n++)
cout << mySignal[n] << ' ';
cout << endl;
return 0;
For completeness, we also include below a MATLAB program that does the same
14
DSP Using MATLAB and Wavelets
thing (assign values to an array, and print them to the screen).
>> mySignal = [ 4.0, 4.6, 5.1, 0.6, 6.0 ];
>> mySignal
mySignal =
4.0000
4.6000
5.1000
0.6000
6.0000
Notice that the MATLAB version is much more compact, without the need to
declare variables first. For simplicity, most of the code in this book is done in
MATLAB.
1.3 Analog Versus Digital
As mentioned earlier, there are two kinds of signals: analog and digital. The word
analog is related to the word analogy; a continuous (“real world”) signal can be
converted to a different form, such as the analog copy seen in Figure 1.3. Here we
see a representation of air pressure sensed by a microphone in time. A cassette tape
recorder connected to the microphone makes a copy of this signal by adjusting the
magnetization on the tape medium as it passes under the read/write head. Thus, we
get a copy of the original signal (air pressure varying with time) as magnetization
along the length of a tape. Of course, we can later read this cassette tape, where the
magnetism of the moving tape affects the electricity passing through the read/write
head, which can be amplified and passed on to speakers that convert the electrical
variations to air pressure variations, reproducing the sound.
An analog signal is one that has continuous values, that is, a measurement can
be taken at any arbitrary time, and for as much precision as desired (or at least as
much as our measuring device allows). Analog signals can be expressed in terms of
functions. A well-understood signal might be expressed exactly with a mathematical
function, or approximately with such a function if the approximation is good enough
for the application. A mathematical function is very compact and easy to work with.
When referring to analog signals, the notation x(t) will be used. The time variable
t is understood to be continuous, that is, it can be any value: t = −1 is valid, as is
t = 1.23456789.
A digital signal, on the other hand, has discrete values. Getting a digital signal
from an analog one is achieved through a process known as sampling, where values
are measured (sampled) at regular intervals, and stored. For a digital signal, the
values are accessed through an index, normally an integer value. To denote a digital
Introduction
15
.
air 
pressure
original
(sound)
time
space (position of tape)
analog copy
magnetization
 .
Figure 1.3: A sound signal with a tape analog.
16
DSP Using MATLAB and Wavelets
signal, x[n] will be used. The variable n is related to t with the equation t = nTs,
where Ts is the sampling time. If you measure the outdoor temperature every hour,
then your sampling time Ts = 1 hour, and you would take a measurement at n = 0
(0 hours, the start time), then again at n = 1 (1 hour), then at n = 2 (2 hours),
then at n = 3 (3 hours), etc. In this way, the signal is quantized in time, meaning
that we have values for the signal only at specific times.
The sampling time does not need to be a whole value, in fact it is quite common
to have signals measured in milliseconds, such as Ts = 0.001 seconds. With this
sampling time, the signal will be measured every nTs seconds: 0 seconds, 0.001
seconds, 0.002 seconds, 0.003 seconds, etc. Notice that n, our index, is still an
integer, having values of 0, 1, 2, 3, and so forth. A signal measured at Ts = 0.001
seconds is still quantized in time. Even though we have measurements at 0.001
seconds and 0.002 seconds, we do not have a measurement at 0.0011 seconds.
Figure 1.4 shows an example of sampling. Here we have a (simulated) continuous
curve shown in time, top plot. Next we have the sampling operation, which is like
multiplying the curve by a set of impulses which are one at intervals of every 0.02
seconds, shown in the middle plot. The bottom plot shows our resulting digital
signal, in terms of sample number. Of course, the sample number directly relates
to time (in this example), but we have to remember the time between intervals for
this to have meaning. In this text, we will start the sample number at zero, just
like we would index an array in C/C++. However, we will add one to the index in
MATLAB code, since MATLAB indexes arrays starting at 1.
Suppose the digital signal x has values x[1] = 2, and x[2] = 4. Can we conclude
that x[1.5] = 3? This is a problem, because there is no value for x[n] when n = 1.5.
Any interpolation done on this signal must be done very carefully! While x[1.5] = 3
may be a good guess, we cannot conclude that it is correct (at the very least, we
need more information). We simply do not have a measurement taken at that time.
Digital signals are quantized in amplitude as well. When a signal is sampled, we
store the values in memory. Each memory location has a finite amount of precision.
If the number to be stored is too big, or too small, to fit in the memory location,
then a truncated value will be stored instead. As an analogy, consider a gas pump.
It may display a total of 5 digits for the cost of the gasoline; 3 digits for the dollar
amount and 2 digits for the cents. If you had a huge truck, and pumped in $999.99
worth of gas, then pumped in a little bit more (say 2 cents worth), the gas pump will
likely read $000.01 since the cost is too large a number for it to display. Similarly,
if you went to a gas station and pumped in a fraction of a cent’s worth of gas, the
cost would be too small a number to display on the pump, and it would likely read
$000.00. Like the display of a gas pump, the memory of a digital device has a finite
amount of precision. When a value is stored in memory, it must not be too large
Introduction
17
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0.18
0.2
−2
−1
0
1
Example continuous signal
0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0.18
0.2
0
0.5
1
Sampling operation
0
1
2
3
4
5
6
7
8
9
10
−2
−1
0
1
2
Resulting digital signal
Figure 1.4: Sampling a continuous signal.
18
DSP Using MATLAB and Wavelets
nor too small, or the amount that is actually stored will be cut to fit the memory’s
capacity. This is what is meant by quantization in amplitude. Incidentally, since
“too small” in this context means a fractional amount below the storage ability of
the memory, a number could be both too large and “too small” at the same time.
For example, the gas pump would show $000.00 even if $1000.004 were the actual
cost of the gas.
Consider the following C/C++ code, which illustrates a point about precision.
unsigned int i;
i=1;
while (i != 0) {
/* some other code */
i++;
}
This code appears to be an infinite loop, since variable i starts out greater than 0 and
always increases. If we run this code, however, it will end. To see why this happens,
we will look at internal representation of the variable i. Suppose our computer has
a word size of only 3 bits, and that the compiler uses this word size for integers.
This is not realistic, but the idea holds even if we had 32, 64, or 128 bits for the
word size instead. Variable i starts out at the decimal value 1, which, of course, is
001 in binary. As it increases, it becomes 010, then 011, 100, 101, 110, and finally
111. Adding 1 results in a 1000 in binary, but since we only use 3 bits for integers
the leading 1 is an overflow, and the 000 will be stored in i. Thus, the condition
to terminate the while loop is met. Notice that this would also be the case if we
remove the unsigned keywords, since the above values for i would be the same, only
our interpretation of the decimal number that it represents will change.
Figure 1.5 shows how a signal can appear in three different ways; as an analog
signal, a digital signal, and as an analog signal based upon the digital version. That
is, we may have a “real world” signal that we digitize and process in a computer,
then convert back to the “real world.” At the top, this figure shows how a simulated
analog signal may appear if we were to view it on an oscilloscope. If we had an
abstract representation for it, we would call this analog signal x(t). In the middle
graph of this figure, we see how the signal would appear if we measured it once every
10 milliseconds. The digital signal consists of a collection of points, indexed by the
sample number n. Thus, we could denote this signal as x[n]. At the bottom of this
figure, we see a reconstructed version of the signal, based on the digital points. We
need a new label for this signal, perhaps x′(t), x̃(t), or x̂(t), but these symbols have
other meanings in different contexts. To avoid confusion, we can give it a new name
altogether, such as x2(t), xreconstructed(t) or new x(t). While the reconstructed
Introduction
19
signal has roughly the same shape as the original (analog) signal, differences are
noticeable. Here we simply connected the points, though other (better) ways do
exist to reconstruct the signal.
0
0.02
0.04
0.06
0.08
0.1
0.12
−2
−1
0
1
2
analog signaltime (seconds)
0
2
4
6
8
10
12
−2
−1
0
1
2
digital signalsample number
0
0.02
0.04
0.06
0.08
0.1
0.12
−2
−1
0
1
2
reconstructed signaltime (seconds)
Figure 1.5: Three ways of viewing a signal.
It is possible to convert from analog to digital, and digital to analog. Analog
signal processing uses electronics parts such as resistors, capacitors, operational
amplifiers, etc. Analog signal processing is cheaper to implement, because these
parts tend to be inexpensive. In DSP, we use multipliers, adders, and delay elements
(registers) to process the signal. DSP is more flexible. For example, error detection
and error correction can easily be implemented in DSP.
1.4 What Is a System?
A system is something that performs an operation (or a transform) on a signal.
A system may be a physical device (hardware), or software. Figure 1.6 shows an
abstraction of a system, simply a “black box” that produces output y[n] based upon
20
DSP Using MATLAB and Wavelets
input x[n]. A simple example system is an incrementer that adds 1 to each value of
the signal. That is, suppose the signal x[n] is {1, 2, 5, 3}. If y[n] = x[n] + 1, then
y[n] would be {2, 3, 6, 4}.
y[n]
x[n]
System
Figure 1.6: An example system.
Examples of systems include forward/inverse Fourier transformers (Chapter 6,
“The Fourier Transform”), and filters such as the Finite Impulse Response filter and
the Infinite Impulse Response filter (see Chapter 3, “Filters”).
1.5 What Is a Transform?
A transform is the operation that a system performs. Therefore, systems and trans-
forms are tightly linked. A transform can have an inverse, which restores the original
values.
For example, there is a cliché that an optimist sees a glass of water as half full
while a pessimist sees it as half empty. The amount of water in the glass does not
change, only the way it’s represented. Suppose a pessimist records the amount of
water in 3 glasses, each of 1 liter capacity, as 50% empty, 75% empty, and 63%
empty. How much water is there in total? To answer this question, it would be
easy to add the amount of water in each glass, but we do not have this information.
Instead, we have measures of how much capacity is left in each. If we put the data
through the “optimist’s transform,” y[n] = 1 − x[n], the data becomes 50%, 25%,
and 37%. Adding these number together results in 112% of a liter, or 1.12 liters of
water.
A transform can be thought of a different way of representing the same informa-
tion.
How full are the glasses of water in Figure 1.7? An optimist would say that they
are half full, one-quarter full, and one-tenth full, respectively.
A pessimist would say the water is 1−(optimist) empty, e.g., .5 empty, .75 empty,
and .90 empty. So we can convert back and forth between the two ways of specifying
the same information. Sometimes one way of looking at the information is better
than another. For example, if you want to know the total amount in all three glasses,
the optimist’s way of seeing things would be easier to work with. However, if you
needed how much more to pour in each glass to fill them, the pessimist’s way of
measuring the amounts is better.
Introduction
21
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
	
































































































































































Figure 1.7: Three glasses of water.
Suppose that we have a simple incrementer system that adds 1 to every input.
We will use x1[n] to represent the input signal, and y1[n] to represent the output. We
can write the output in terms of the input, i.e., y1[n] = x1[n]+1. If x1 is the sequence
{1, 3, 7}, the output y1 would be the sequence {2, 4, 8}. Now suppose we want to
get the original values of x1 back. To do this, we need to “undo” the transform,
which we can do with the following “inverse” transform. Since we have given a
name to the reverse transform, we will call the first transform (y1[n] = x1[n] + 1)
the “forward” transform.
The inverse transform for the incrementer system would be a decrementer, that
simply subtracts 1 from each value, i.e., y2[n] = x2[n] − 1. If we hooked these two
systems up in series, we would have the output of the first system as the input to the
second, effectively assigning x2[n] = y1[n]. So y2[n] = x2[n] − 1, y2[n] = y1[n] − 1,
or y2[n] = (x1[n] + 1) − 1. It is easy to see that the two systems cancel each other
out, and the final output (y2[n]) is what we started with originally (x1[n]).
Typically, we do not study systems that add or subtract constants.
A natural question is, “Why would we do this?” One answer is so that we can
analyze the transformed signal, or compress it. This does not follow from such a
simple example, but more complicated transforms (such as the discrete cosine trans-
form) have been used effectively to allow a compression program to automatically
alter the signal to store it in a compact manner.
Sometimes transforms are performed because things are easier to do in the trans-
formed domain. For example, suppose you wanted to find the number of years
between the movies Star Wars (copyright MCMLXXVII) and Star Wars: Episode
II—Attack of the Clones (copyright MMII). Subtracting one number from another
is not a problem, but Roman numerals are difficult to work with [7]! The algorithm
we use of subtracting the right-most column digits first, writing the difference and
borrowing from the left, does not work. Instead, the easiest thing to do is convert
the Roman numerals to decimal numbers, perform the subtraction, and convert the
result back to Roman numerals. In this case, the transform would convert MCM-
LXXVII to 1977, and MMII to 2002. Subtracting 1977 from 2002 gives us 25, which
is the inverse-transformed to XXV. Performing the transform and later the inverse-
22
DSP Using MATLAB and Wavelets
transform may seem like a waste of effort, but not if it greatly simplifies the steps
in between.
If you have had a course in differential equations, you may have seen the Laplace
transform used for just this purpose: to provide an easy way to solve these equations.
After the equation is set up, you apply the Laplace transform (often by looking it
up in a table), then the transformed equation can be manipulated and solved using
algebra, and finally the result is inverse-transformed back to the time-domain (again,
often by looking at a table). The Laplace transform also has a use in determining
the stability of systems, and this is something we will return to in a later chapter.
1.6 Why Do We Study Sinusoids?
We look at sinusoidal functions (sine and cosine) because they are interesting func-
tions that often appear in the “analog” world. Examining a single cosine function is
easy to do, and what applies to one cosine function also applies to a signal composed
of several sinusoids. Since sinusoidal functions occur frequently, it is nice to be able
to use a few pieces of data to represent a sinusoid. Almost every sinusoid in this
text will be of the following format:
amplitude × cos(2π× frequency × t + phase).
That is all! When the amplitude, frequency, and phase are known, we have all the
information we need to find the value of this sinusoid for any value of time (t). This
is a very compact way of representing a signal. If the signal is more complex, say
it is composed of two sinusoids, we just need the amplitude, frequency, and phase
information of the two sinusoids. In this way, we can represent the signal, then later
remake it from this information.
Think of how the earth rotates around the sun in a year’s time. The earth gets
the same amount of sunlight every day, but the earth’s tilt means that one hemi-
sphere or the other receives more sun. Starting at the December solstice (around
December 22), the earth’s Northern hemisphere receives the least amount of sunlight
[8]. As the earth rotates around the sun, this hemisphere gets an increasing amount
of sunlight (spring time), until it reaches the maximum (summer), then the days get
shorter (fall), and it finally returns to the low point again in winter. The seasons
vary with time. In two dimensions, one could diagram this with the sun in the
center, and the earth following a circular path around it (approximately). Another
way of viewing this information is to graph the distance with time as the x-axis.
For the y-axis, consider the horizontal distance between the two. This graph looks
like a sine wave. The seasons change with time, but they always repeat. Therefore,
this signal is periodic.
Introduction
23
Examples of sinusoidal signals are all around us. The rotation of the moon
around the earth, the temperature of a city over a day’s time, and the acoustic
waves made by human speech are all examples of sinusoidal-like signals. Human
speech is interesting, because it may be a single sinusoid (such as a vowel sound),
or a composite signal made up of several sinusoids.
Composite signals are important in theory as well as practice. Any function that
is not truly random can be approximated with a sum of sinusoids. In other words, a
signal can be broken down into a combination of simpler signals. As a mathematical
abstraction, the composite signal
x(t) = cos(2π100t) + cos(2π200t)
can easily be plotted in MATLAB. The example below will work, though it is inef-
ficient.
Ts = 0.0002;
% sampling time, Ts, is .2 ms
for n = 1:100
% we will have 100 samples
x(n) = cos(2*pi*100*n*Ts) + cos(2*pi*200*n*Ts);
end
plot(x);
A better way is to eliminate the for loop, as shown below.
Ts = 0.0002; % sampling time, Ts, is .2 ms
n = 1:100;
% we will have 100 samples
x = cos(2*pi*100*n*Ts) + cos(2*pi*200*n*Ts);
plot(n, x, 'b');
The code above adds the two cosine signals together, point by point, and stores
the results in a vector x [6] [9] [10]. We can use n in place of 1:length(x) for this
particular example, if we want. Also, using x would work here just as well as x(n).
Notice that the MATLAB example above does not actually plot x(t). It plots
x(t = nTs), where n is an integer 1..1000. Ts is the sampling time, which is an
important value that we will look into more closely in a later chapter. We could
simulate x(t), however, it would still be a digital signal since variable t has a limited
set of values. Example code appears below.
% simulated x(t)
t = 0:0.0002:0.02;
x = cos(2*pi*100*t) + cos(2*pi*200*t);
plot(1:length(x), x, 'b');
24
DSP Using MATLAB and Wavelets
The 1:length(x) parameter specifies an array that has the same number of
elements as x does, which can be seen along the x-axis. A final example of the plot
command follows:
plot(t, x, 'b');
Notice how the x-axis now reflects time instead of sample number.
When several sinusoid components are in the same signal, and each of the fre-
quencies is an integer multiple of a common frequency, we say that the frequencies
are harmonically related.
A “truly random” function cannot be represented well with a sum of sinusoids.
Such a function would not repeat itself, but would go on forever as an unpredictable
value for any time value. This is similar to an irrational number. Most values are
rational, that is, can be represented as a division of one integer number by another,
for example, the quantity 1.5 is a rational number since it is the same as 3/2. It
could also be represented as 15/10, if desired, just as 1.4 could be represented as
14/10, and 1.3 could be represented as 13/10. In fact, any number that has only
1 digit past the decimal point could be written as an integer over 10. In a similar
line of reasoning, any number with 2 digits past the decimal could be written as an
integer over 100, e.g., 1.23 = 123/100. We can keep going with this, even with an
extreme case like 1.23456789 = 123456789/10000000. Since our denominator (the
number on the bottom) can be as large as we like, it may seem like EVERY number
can be represented in such a way, and is therefore rational. But there are a few
numbers that are irrational, such as π. You can write as many digits for π as you
like, and someone could come along after you and add another digit. It does not
repeat, as far as we know, even when computed out to millions of digits. Truly
random functions are like this—they do exist, but we can deal with them the way
we deal with π, by approximation. One could represent π as the rational number
3.14159. This is not TRULY π, but it is close enough for most applications. In
a similar fashion, one could find an approximation for a truly random signal that
would be judged “good enough” for the application.
1.7 Sinusoids and Frequency Plots
Suppose you had the following sequence of numbers, with little information about
them: {59, 65, 68, 70, 61, 61, 63, 66, 66, 67, 70, 71, 71, 70, 70, 70}. You could
calculate the average, the variance, and what the minimum and maximum values
are. If you had to guess what the signal is likely to be, measured at an arbitrary
time, 66 seems to be a good guess. Now suppose you get the information that these
Introduction
25
values correspond to temperature readings in degrees Fahrenheit. These numbers
could be used to decide how to dress for the climate, and you might conclude that
long-sleeved garments would be warm enough, perhaps with a sweater. Now suppose
that these readings were taken in Atlanta during the month of August. How can
this be? Atlanta is known to have a hot climate, especially in August. There is one
key piece of information that is missing: that these temperatures were taken every
day at 2 a.m. Now it all makes sense, that the temperatures were read at a time
when they are likely to be their lowest.
Taking samples too infrequently leads to poor conclusions. Here one might
assume that the temperatures fluctuate by only a few degrees. Or one could as-
sume that the temperatures above are indicative of the climate, instead of being
an extreme. We expect that temperature will rise in the afternoon and fall in the
nighttime, not unlike a sine function. But this example shows how we can miss the
overall rise and fall trend by sampling at regular intervals that happen to coincide
with the low points. If we were to sample a bit more frequently, say every 23 hours
starting with an initial read at 2 a.m., this new set of readings would not be better.
They would make it appear as if the temperature gradually gets warmer (as the
sampling time becomes earlier: first 2 a.m., then 1 a.m., then midnight, then 11
p.m., etc., eventually reaching the afternoon), then gradually gets cooler (when the
samples get earlier and earlier in the morning). A better sampling rate would be
every 12 hours or less, since we would get a much better idea of how the temper-
ature fluctuates. At worst, we would read temperatures right in the middle of the
extremes, but we would still have a good idea of the average.
We use both time and frequency in language to view the same information in
different ways. For example, you may ask “how long is it before I have to change
the oil in my car?” in which case the answer would be in terms of time (months).
A related question, “how many times a year do I need to change the oil in my car?”
results in a frequency-based answer (oil changes per year). Time and frequency have
an inverse relationship in language as well as in DSP.
We are used to thinking about frequency in some contexts. When you move to
a new town, you are likely to set up your stereo and scan the channels. Moving the
dial as far as it will go to the left, you slowly turn it until you hear something. If
you like what you hear, you make note of the dial’s location, and move it to the
right. When finished, you have essentially made a graph of the spectrum: where
you found nonzero energy for the dial’s location, you marked the frequency.
You may wonder why radio stations often have numbers past the decimal point,
such as 88.5 instead of just 88 or 89. These numbers are in terms of MHz, or millions
of Hertz, and the extra digit helps you tune in a good signal. After you graduate
and get a great job, and see that you have $1.3M in your bank account, think back
26
DSP Using MATLAB and Wavelets
to this example, and how the .3 is significant to you! You may have noticed how FM
radio stations end in odd numbers, at least in the United States. The U.S. allocated
the FM spectrum in 200 kHz chunks, starting on an odd boundary. So you might
have a local 88.5 MHz station, but not one at 88.6 MHz.
Now imagine instead that you have a measure of all the energy of the radio
spectrum for an instant in time. How could you use this to tell which radio station
to listen to? You cannot use it for this purpose. Though you are used to signals
in terms of time, there are some cases in which the frequency information makes
things much easier.
1.8 Summations
In Digital Signal Processing, you are likely to run into formulas containing summa-
tions and other seemingly complicated things. Summations are nothing to be afraid
of. For example, the following formula and the C/C++ code below it accomplish
the same task. This example summation is just a compact way of expressing an
idea: for each number in x, multiply by two and subtract one, add all the results
together, and store the final answer in the variable s. Let x = {1.4, 5.7, 3.0, 9.2,
6.8}, and let N be the number of values in x.
s =
N−1
∑
i=0
2x[i] − 1
Here is part of a program in C/C++ to do the same thing:
int i;
float x[] = {1.4, 5.7, 3.0, 9.2, 6.8};
int N=5; // the number of values in x[]
float s = 0;
for (i=0; i<N; i++)
s = s + 2*x[i] - 1;
If you find the mathematical notation confusing, then practice a few examples of
converting sums to programs. After a while, you should be comfortable with both
ways.
In MATLAB, we could do the following:
x = [1.4, 5.7, 3.0, 9.2, 6.8];
Introduction
27
s = 0;
for i=1:length(x)
s = s + 2*x(i) - 1;
end
Even more compactly, we could accomplish the summation in MATLAB with:
x = [1.4, 5.7, 3.0, 9.2, 6.8];
s = sum(2*x - 1);
MATLAB will be covered in more detail in Chapter 2.
1.9 Summary
This chapter covers the concepts of numbers, from counting to a number line, to neg-
ative numbers, infinity and negative infinity, fractional numbers, irrational numbers,
and finally complex numbers (
√
−1 being the y-axis). Topics also covered include
polar coordinates, complex exponentials, continuous versus discrete, digital versus
analog, signals, systems, transforms, sinusoids, and composite signals. Coverage of
these topics will be expanded in the following chapters.
Signals are the data that we work with, and digital signals can be thought of as
an array (or matrix). Online systems would have data in streams, but we typically
do not discuss these to keep things simple. Transforms can be used for analysis, and
are often used in compression applications.
1.10 Review Questions
1. What is a signal?
2. What is a system?
3. What is a transform?
4. What does x[n] imply?
5. What does x(t) imply?
28
DSP Using MATLAB and Wavelets
6. What is the decimal value of the binary fixed-point number 0101.01? (show
work)
7. Represent 75.625 as a fixed-point binary number.
8. Represent 36.1 as a fixed-point binary number.
9. Convert the Cartesian point (2, 5) to polar coordinates.
10. Convert the polar coordinates 4, 25 degrees to Cartesian coordinates.
11. Explain what we mean by quantization in time, and quantization in amplitude.
12. Why are sinusoids useful in digital signals?
13. Why is there a problem using θ = arctan(b/a) for conversion to polar coordi-
nates?
14. What are the differences between analog and digital?
Chapter 2
MATLAB
MATLAB is a programming environment which has gained popularity due to its
ease of use. The purpose of this chapter is not to teach MATLAB programming,
but it is intended to provide a quick reference to MATLAB. It is easy to understand
if you already know a programming language; if you do not, the following examples
should help. Whether you are an experienced programmer or a beginner, the best
thing to do is try the examples on your computer. MATLAB provides an interactive
environment that makes this easy: you type a command, and it executes it. If there
is a problem, it will let you know, and you can try again. If a command works well,
you can copy and paste it into a file, creating a complex program one step at a time.
In fact, we will consider any group of commands that are stored in a file to be a
program. The only distinction between a program and a function is the “function”
keyword used in the first (noncomment) line in a file. This line specifies the inputs
and outputs, as well as the function name. It is a good idea to give the function
the same name as the file that it is stored in. The filename should also have “.m”
as the extension. An example function will follow shortly.
MATLAB provides an editor, which will color text differently based on its func-
tionality. Under this editor, comments appear in green, strings appear in red, and
keywords (such as if, end, while) appear in blue. This feature helps the program-
mer spot errors with a command, such as where an unterminated string begins.
This editor is not required, however, because any text editor will work. Here is one
frustrating mistake that programmers occasionally make: if a program or function
is changed in the editor, but not saved to disk, then the old version of it will be used
when it is executed. Make sure to save your work before trying to run it!
29
30
DSP Using MATLAB and Wavelets
2.1 Working with Variables
A common programming mistake is to use a variable without declaring it first. This
is not a problem in MATLAB. In fact, MATLAB is very flexible about variables
and data types. If an undefined variable is assigned a value, for example:
a = 1;
MATLAB will assume that the variable a should be created, and that the value 1
should be stored in it. If the variable already exists, then it will have the number 1
replace its previous value.
This works with scalars, vectors, and matrices. The numbers stored in variables
can be integers (e.g., 3), floating-points (e.g., 3.1), or even complex (e.g., 3.1+5.2j).
a = 3;
b = 3.1;
c = 3.1+5.2j;
A semicolon “;” is often put at the end of a command. This has the effect of
suppressing the output of the resulting values. Leaving off the semicolon does not
affect the program, just what is printed to the screen, though printing to the screen
can make a program run slowly. If a variable name is typed on a line by itself,
then the value(s) stored for that variable will be printed. This also works inside a
function or a program. A variable on a line by itself (without a semicolon) will be
displayed.
>> a
a =
3
>> b
b =
3.1000
>> c
MATLAB
31
c =
3.1000 + 5.2000i
MATLAB shows the complex value stored in c with i representing the complex
part (the square root of −1). Notice, though, that it also understands j, as the
assignment statement for c shows. Both i and j can be used as variables, just as
they are in other languages. It is advisable to use k instead, or some other variable
name, to avoid confusion.
It is a good idea to use uppercase variable names for matrices, and lowercase
names for scalars. This is not a requirement; MATLAB does not care about this.
But it is case sensitive; that is, the variable x is considered to be different from the
variable X.
One final note about functions is that the variables used in a function are not
accessible outside of that function. While this helps the programmer maintain local
variables, it also means that, if the function produces an error, the programmer
will not be able to simply type a variable’s name at the prompt to find its value.
A quick fix for this is to comment out the function’s declaration line, then call
it again. Alternatively, MATLAB provides breakpoints and other programming
interface features that you would expect from a software development kit.
2.2 Getting Help and Writing Comments
The help command is very useful. Often, a programmer will remember the function
he or she wants to use, but first want to verify the name, and then figure out what
the parameters are, and the order of the parameters. The help command is useful
for this purpose. Just type help at the prompt, followed by the function name, for
example, help conv gives information on the convolution function (conv).
When creating a function or program, it is important to document the code.
That is, make sure to include enough comments so that someone else can figure
out what the code does. This “someone else” could even be yourself, if you make
a program and enough time passes before you need to use it again. Comments in
MATLAB start with a percent sign “%,” and the comment lasts until the end of
the line, just like the C++ style comment (of two forward slashes). A comment can
appear on the same line as a command, as long as the command is first (otherwise,
it would be treated as part of the comment).
It is a good idea to put several lines of comments at the top of your programs and
functions. Include information like the order of the inputs, the order of the outputs,
32
DSP Using MATLAB and Wavelets
and a brief description of what the code does. The nice thing about MATLAB is
that the comments at the top of your file will be returned by the help function.
Try this: Type help xyz at the MATLAB prompt. It should tell you that no
such function exists.
>> help xyz
xyz.m not found.
Type the function below, and save it as “xyz.m”:
%
% This is my function. It returns the input as the output.
%
%
Usage:
a = xyz(b)
%
function [a] = xyz(b)
a=b;
Now try help xyz again. You should see the first few lines printed to the screen.
Be sure to put the percent sign in the left-most column if you want these comments
to be found by the help facility.
>> help xyz
This is my function. It returns the input as the output.
Usage:
a = xyz(b)
2.3 MATLAB Programming Basics
This section covers the main things that you need to know to work with MATLAB.
The command clc clears the command window, and is similar to the clear com-
mand in the Unix environment. The command clear, by the way, tells MATLAB