A Primer on WAVELETS and Their Scientific Applications SECOND EDITION
© 2008 by Taylor & Francis Group, LLC
C7451_FM.i...

Author:
James S. Walker

This content was uploaded by our users and we assume good faith they have the permission to share this book. If you own the copyright to this book and it is wrongfully on our website, we offer a simple DMCA procedure to remove your content from our site. Start by pressing the button below!

A Primer on WAVELETS and Their Scientific Applications SECOND EDITION

© 2008 by Taylor & Francis Group, LLC

C7451_FM.indd 1

1/3/08 11:28:56 AM

A Primer on WAVELETS and Their Scientific Applications SECOND EDITION

James S. Walker University of Wisconsin Eau Claire, Wisconsin, U.S.A.

© 2008 by Taylor & Francis Group, LLC

C7451_FM.indd 3

1/3/08 11:28:56 AM

Chapman & Hall/CRC Taylor & Francis Group 6000 Broken Sound Parkway NW, Suite 300 Boca Raton, FL 33487-2742 © 2008 by Taylor & Francis Group, LLC Chapman & Hall/CRC is an imprint of Taylor & Francis Group, an Informa business No claim to original U.S. Government works Printed in the United States of America on acid-free paper 10 9 8 7 6 5 4 3 2 1 International Standard Book Number-13: 978-1-58488-745-4 (Softcover) This book contains information obtained from authentic and highly regarded sources Reasonable efforts have been made to publish reliable data and information, but the author and publisher cannot assume responsibility for the validity of all materials or the consequences of their use. The Authors and Publishers have attempted to trace the copyright holders of all material reproduced in this publication and apologize to copyright holders if permission to publish in this form has not been obtained. If any copyright material has not been acknowledged please write and let us know so we may rectify in any future reprint Except as permitted under U.S. Copyright Law, no part of this book may be reprinted, reproduced, transmitted, or utilized in any form by any electronic, mechanical, or other means, now known or hereafter invented, including photocopying, microfilming, and recording, or in any information storage or retrieval system, without written permission from the publishers. For permission to photocopy or use material electronically from this work, please access www. copyright.com (http://www.copyright.com/) or contact the Copyright Clearance Center, Inc. (CCC) 222 Rosewood Drive, Danvers, MA 01923, 978-750-8400. CCC is a not-for-profit organization that provides licenses and registration for a variety of users. For organizations that have been granted a photocopy license by the CCC, a separate system of payment has been arranged. Trademark Notice: Product or corporate names may be trademarks or registered trademarks, and are used only for identification and explanation without intent to infringe. Library of Congress Cataloging-in-Publication Data Walker, James S. A primer on wavelets and their scientific applications / James S. Walker. -- 2nd ed. p. cm. Includes bibliographical references and index. ISBN 978-1-58488-745-4 (alk. paper) 1. Wavelets (Mathematics) I. Title. QA403.3.W33 2008 515’.2433--dc22

2007048858

Visit the Taylor & Francis Web site at http://www.taylorandfrancis.com and the CRC Press Web site at http://www.crcpress.com © 2008 by Taylor & Francis Group, LLC

C7451_FM.indd 4

1/3/08 11:28:56 AM

i

i

i

i

To Berlina and Damian

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

About the author James S. Walker received his doctorate from the University of Illinois, Chicago, in 1982. He has been a professor of Mathematics at the University of Wisconsin-Eau Claire since 1982. He has published papers on Fourier analysis, wavelet analysis, complex variables, and logic, and is the author of four books on Fourier analysis, FFTs, and wavelet analysis. He has been married to Angela Huang since 1993, and they are the proud parents of two children, Berlina and Damian, to whom this book is dedicated.

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

Preface

“Wavelet theory” is the result of a multidisciplinary eﬀort that brought together mathematicians, physicists and engineers...this connection has created a ﬂow of ideas that goes well beyond the construction of new bases or transforms. —St´ephane Mallat1

The past two decades have witnessed an explosion of activity in wavelet analysis. Thousands of research papers have been published on the theory and applications of wavelets. Wavelets provide a powerful and remarkably ﬂexible set of tools for handling fundamental problems in science and engineering. With the second edition of this primer, I aim to provide the most up-to-date introduction to this exciting ﬁeld. What is new in this second edition? I have added a lot of new material to this second edition. The most important additions are the following: • Exercises and additional examples: At the end of the chapters, I provide a total of 104 worked examples, and 222 exercises (with solutions provided for representative problems). These examples and exercises constitute a book-in-itself of review material, and should make this primer a valuable text for courses on wavelets. • Biorthogonal wavelets: I have added two sections that describe the most important examples of these wavelets, and I describe their application to image compression. 1 Mallat’s

quote is from [4] listed in the Notes and references section for Chapter 3.

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

• Image compression: I have included a mini-course on image compression, with descriptions of two state-of-the art algorithms (ASWDR and JPEG 2000). This discussion includes a tutorial on arithmetic compression, an important topic that is omitted from most (if not all) introductions to image compression. • Image denoising: I have included a lot more material on image denoising, including a description of a high-performance wavelet denoiser. In addition to the standard topic of removing Gaussian random noise, which most books limit themselves to, I also describe a technique for removing isolated, randomly positioned, clutter from images. • Gabor transforms and time-frequency analysis: I provide a concise introduction to the fundamentals of time-frequency analysis and show how it is applied to musical theory, musical synthesis, and audio denoising. An important new mathematical concept in musical theory, known as the multiresolution principle, is described here for the ﬁrst time in book form. I also provide an extensive list of references to the ﬁeld of timefrequency analysis, a topic that is given short shrift in most introductions to wavelets. • Project ideas: I provide a list, with brief descriptions, of a set of research projects on wavelets and time-frequency analysis. See Appendix A. • Book website: To keep the book as current as possible, and to provide a wealth of additional online material such as software, sound and image ﬁles, and links to other web resources on wavelets, I have created a website for the book. I describe this website in more detail later in the preface. • Enhanced list of references: I have more than doubled the list of references for the book; it now totals about 150 items. At the website for the book I provide an online version of this bibliography with a great many links for free downloading of important papers on wavelets and time-frequency analysis. What problems can be solved with wavelets? For an idea of the wide range of problems that can be solved with wavelets, here is a list of some of the problems that I shall discuss: • Audio denoising: Long distance telephone messages, for instance, often contain signiﬁcant amounts of noise. How do we remove this noise? • Signal compression: The eﬃcient transmission of large amounts of data, over the Internet for example, requires some kind of compression. How do we compress this data without losing signiﬁcant information? © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

• Object detection: What methods can we use to pick out a small image, say of an aircraft, from the midst of a larger more complicated image? • Image denoising: Images formed by microscopes and other devices are often contaminated by large amounts of unwanted clutter (referred to as noise). Can this noise be removed in order to clarify the image? • Image enhancement: When an optical microscope image is recorded, it often suﬀers from blurring. How can the appearance of the objects in these images be sharpened? • Image recognition: How do humans recognize faces? Can we teach machines to do it? • Diagnosing heart trouble: Is there a way to detect abnormal heartbeats, hidden within a complicated electrocardiogram? • Speech recognition: What factors distinguish consonants from vowels? How do humans recognize diﬀerent voices? • Musical analysis. How can mathematics help us better understand the nature of musical composition? All of these problems can be tackled using wavelets. I will show how during the course of this book. Why do we need another wavelet book? The goal of this primer is to guide the reader through the main ideas of wavelet analysis in order to facilitate a knowledgeable reading of the present research literature, especially in the applied ﬁelds of audio and image processing and biomedicine. Although there are many excellent books on the theory of wavelets, these books are focused on the construction of wavelets and their mathematical properties. Furthermore, they are all written at a graduate school level of mathematical and/or engineering expertise. There remains a real need for a basic introduction, a primer, which uses only elementary algebra and a smidgen of calculus to explain the underlying ideas behind wavelet analysis, and devotes the majority of its pages to explaining how these underlying ideas can be applied to solve signiﬁcant problems in audio and image processing and in biology and medicine. To keep the mathematics as elementary as possible, I focus on the discrete theory. It is in the continuous theory of wavelet analysis where the most diﬃcult mathematics lies; yet when this continuous theory is applied it is almost always converted into the discrete approach that I describe in this primer. Focusing on the discrete case will allow us to concentrate on the applications of wavelet analysis while at the same time keeping the mathematics under control. On the rare occasions when we need to use more advanced mathematics, I shall mark these discussions oﬀ from the main text by putting them © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

into subsections that are marked by asterisks in their titles. An eﬀort has been made to ensure that subsequent discussions do not rely on this more advanced material. Chapter summaries Chapter 1 summarizes our subject. Using examples from imaging and audio, I brieﬂy explain wavelet analysis in non-mathematical terms. In Chapter 2 I introduce the simplest wavelets, the Haar wavelets. I also introduce many of the basic concepts—wavelet transforms, energy conservation and compaction, multiresolution analysis, compression and denoising—that will be used in the remainder of the book. For this reason, I devote more pages to the theory of Haar wavelets than perhaps they deserve alone; keep in mind that this material will be ampliﬁed and generalized throughout the remainder of the book. In Chapter 3 I describe the Daubechies wavelets, which have played a key role in the explosion of activity in wavelet analysis. After a brief introduction to their mathematical properties, I describe several applications of these wavelets. First, I explain in detail how they can be used to compress audio signals—this application is vital to the ﬁelds of telephony and telecommunications. Second, I describe how the method of thresholding provides a powerful technique for removing random noise (static) from audio signals. Removing random noise is a fundamental necessity when dealing with all kinds of data in science and engineering. The threshold method, which is analogous to how our nervous system responds only to inputs above certain thresholds, provides a nearly optimal method for removing random noise. Besides random noise, Daubechies wavelets can also be used to remove isolated “pop-noise” from audio. The chapter concludes with an introduction to biorthogonal wavelets, which have played an important part in wavelet image compression algorithms such as JPEG 2000. Wavelet analysis can also be applied to images. In Chapter 4 I provide a mini-course on wavelet-based image compression. I describe two state-of-theart image compression algorithms, including the new JPEG 2000 standard. I also discuss wavelet-based image denoising, including examples motivated by optical microscopy and scanning tunnelling microscopy. Chapter 4 concludes with some examples from image processing. I examine edge detection, and the sharpening of blurred images, and an example from computer vision where wavelet methods can be used to enormously increase the speed of identiﬁcation of an image. Chapter 5 relates wavelet analysis to frequency analysis. Frequency analysis, also known as Fourier analysis, has long been one of the cornerstones of the mathematics of science and engineering. I shall brieﬂy describe how wavelets are characterized in terms of their eﬀects on the frequency content of signals. One application that I discuss is object identiﬁcation—locating a small object within a complicated scene—where wavelet analysis in concert © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

with Fourier analysis provides a powerful approach. The chapter concludes with a short introduction to the theory and application of Gabor transforms, an important relative of wavelet transforms. In recent years Gabor transforms have been increasingly used in sound processing and denoising, where they may be even more useful than wavelet transforms. In the ﬁnal chapter I deal with some extensions which reach beyond the fundamentals of wavelets. I describe a generalization of wavelet transforms known as wavelet packet transforms. I apply these wavelet packet transforms to compression of audio signals, images, and ﬁngerprints. Then I turn to the subject of continuous wavelet transforms, as they are implemented in a discrete form on a computer. Continuous wavelet transforms are widely used in seismology and have also been used very eﬀectively for analyzing speech and electrocardiograms. We shall also see how they apply to analyzing musical rhythm. Unfortunately, for reasons of space, I could not examine several important aspects of wavelet theory, such as local cosine series, best-basis algorithms, or fast matrix multiplication. At the end of the ﬁnal chapter, in its Notes and references section, I provide a guide to the literature for exploring omitted topics. Computer software Without question the best way to learn about applications of wavelets is to experiment with making such applications. This experimentation is typically done on a computer. In order to simplify this computer experimentation, I have created free software, called FAWAV. Further details about FAWAV can be found in Appendix C; suﬃce it to say for now that it is designed to allow the reader to duplicate all of the applications described in this primer and to experiment with other ideas. In this second edition I have added several new features, including image compression and denoising of both images and audio, which should add to the usefulness of FAWAV for digital signal processing. Notes and references In the Notes and references sections that conclude each chapter, I provide the reader with ample references where further information can be found. These references are numbered consecutively within each chapter. At the risk of some duplication, but with the reader’s convenience in mind, I have also compiled all references into a comprehensive Bibliography. Primer website To keep the material in the book as up-to-date as possible, and to facilitate classroom use, there is a website for this book at the following address: http://www.uwec.edu/walkerjs/Primer/ © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

This website contains a lot of supplementary material for the book, including the following: • The latest version of FAWAV and its User Manual. The User Manual is a comprehensive introduction to the use of FAWAV for digital signal processing. • Papers for downloading, including links to all papers cited in the primer. • Sound and image ﬁles for use with FAWAV, including all the sound and image ﬁles referred to in the primer. • A collection of Project ideas for carrying out your own research. This online collection will include an archive of preprints and publications that deal with these projects, and supplement the list from the primer with new topics. • Links to other websites dealing with wavelet and time-frequency analysis, including all the websites listed in the primer. Acknowledgments It is a pleasure to acknowledge several people who have helped me in writing this book. First, I must thank my wife Angela for her constant love and support. I would also like to thank my executive editor, Bob Stern, for his encouragement and help. I appreciate several suggestions from Steve Krantz that helped me put some ﬁnishing touches on the book. Thanks to Thomas Maresca for some enlightening discussion on FAWAV and its use in this primer. Thanks also to the UWEC Oﬃce of University Research and Creative Activities for a grant that gave me some extra time away from my teaching responsibilities for work on this project. Finally, my students in my Digital Signal Processing and Digital Image Processing courses deserve a warm thank you for all of their many suggestions and questions, which helped me to write a better book. I especially would like to thank Xiaowen Cheng, Jarod Hart, Nicholas Nelson, Derek Olson, and Elizabeth Otteson for their contributions. James S. Walker Eau Claire, Wisconsin

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

Contents 1 Overview 1.1 What is wavelet analysis? . . . . . . . . . . . . . . . . . . . . . 1.2 Notes and references . . . . . . . . . . . . . . . . . . . . . . . . 2 Haar wavelets 2.1 The Haar transform . . . . . . . . . . . . . . . 2.1.1 Haar transform, 1-level . . . . . . . . . 2.2 Conservation and compaction of energy . . . . 2.2.1 Conservation of energy . . . . . . . . . . 2.2.2 Haar transform, multiple levels . . . . . 2.2.3 Justiﬁcation of conservation of energy . 2.3 Haar wavelets . . . . . . . . . . . . . . . . . . . 2.4 Multiresolution analysis . . . . . . . . . . . . . 2.4.1 Multiresolution analysis, multiple levels 2.5 Signal compression . . . . . . . . . . . . . . . . 2.5.1 A note on quantization . . . . . . . . . 2.6 Removing noise . . . . . . . . . . . . . . . . . . 2.6.1 RMS Error . . . . . . . . . . . . . . . . 2.7 Notes and references . . . . . . . . . . . . . . . 2.8 Examples and exercises . . . . . . . . . . . . .

1 1 4

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

5 6 7 9 10 11 12 14 16 19 21 26 26 29 30 31

3 Daubechies wavelets 3.1 The Daub4 wavelets . . . . . . . . . . . . . . . . . . 3.1.1 Remarks on small ﬂuctuation values ∗ . . . . 3.2 Conservation and compaction of energy . . . . . . . 3.2.1 Justiﬁcation of conservation of energy ∗ . . . 3.2.2 How wavelet and scaling numbers are found ∗ 3.3 Other Daubechies wavelets . . . . . . . . . . . . . . 3.3.1 Coiﬂets . . . . . . . . . . . . . . . . . . . . . 3.4 Compression of audio signals . . . . . . . . . . . . . 3.4.1 Quantization and the signiﬁcance map . . . . 3.5 Quantization, entropy, and compression . . . . . . . 3.6 Denoising audio signals . . . . . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

41 41 49 50 50 53 54 58 61 62 65 69

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

3.6.1 Choosing a threshold value . . . . . 3.6.2 Removing pop noise and background 3.7 Biorthogonal wavelets . . . . . . . . . . . . 3.7.1 Daub 5/3 system . . . . . . . . . . . 3.7.2 Daub 5/3 inverse . . . . . . . . . . . 3.7.3 MRA for the Daub 5/3 system . . . 3.7.4 Daub 5/3 transform, multiple levels 3.7.5 Daub 5/3 integer-to-integer system . 3.8 The Daub 9/7 system . . . . . . . . . . . . 3.9 Notes and references . . . . . . . . . . . . . 3.10 Examples and exercises . . . . . . . . . . .

. . . . static . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

70 73 75 76 78 78 80 82 83 85 87

4 Two-dimensional wavelets 4.1 Two-dimensional wavelet transforms . . . . . . . . . . . . 4.1.1 Discrete images . . . . . . . . . . . . . . . . . . . . 4.1.2 2D wavelet transforms . . . . . . . . . . . . . . . . 4.1.3 2D wavelets and scaling images . . . . . . . . . . . 4.2 Compression of images—fundamentals . . . . . . . . . . . 4.2.1 Error measures in image processing . . . . . . . . . 4.2.2 Comparing JPEG with wavelet-based compressors 4.3 Fingerprint compression . . . . . . . . . . . . . . . . . . . 4.4 The WDR algorithm . . . . . . . . . . . . . . . . . . . . . 4.4.1 Bit-plane encoding . . . . . . . . . . . . . . . . . . 4.4.2 Diﬀerence reduction . . . . . . . . . . . . . . . . . 4.4.3 Arithmetic compression . . . . . . . . . . . . . . . 4.5 The ASWDR algorithm . . . . . . . . . . . . . . . . . . . 4.5.1 Arithmetic compression . . . . . . . . . . . . . . . 4.5.2 Relation to vision . . . . . . . . . . . . . . . . . . . 4.6 Important image compression features . . . . . . . . . . . 4.6.1 Progressive transmission/reconstruction . . . . . . 4.6.2 Lossless compression . . . . . . . . . . . . . . . . . 4.6.3 Region-of-interest . . . . . . . . . . . . . . . . . . . 4.7 JPEG 2000 image compression . . . . . . . . . . . . . . . 4.7.1 Compressing color images . . . . . . . . . . . . . . 4.8 Denoising images . . . . . . . . . . . . . . . . . . . . . . . 4.8.1 The TAWS algorithm . . . . . . . . . . . . . . . . 4.8.2 Comparison with Wiener denoising . . . . . . . . . 4.8.3 Estimation of noise standard deviation ∗ . . . . . . 4.8.4 Removal of clutter noise . . . . . . . . . . . . . . . 4.9 Some topics in image processing . . . . . . . . . . . . . . 4.9.1 Edge detection . . . . . . . . . . . . . . . . . . . . 4.9.2 Edge enhancement . . . . . . . . . . . . . . . . . . 4.9.3 Image recognition . . . . . . . . . . . . . . . . . . 4.10 Notes and references . . . . . . . . . . . . . . . . . . . . . 4.11 Examples and exercises . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

97 97 98 99 102 104 107 108 110 113 113 116 119 123 125 126 127 127 128 130 130 132 133 133 134 136 137 139 139 140 141 144 147

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

5 Frequency analysis 5.1 Discrete Fourier analysis . . . . . . . . . . . . . . . . . 5.1.1 Frequency content of wavelets . . . . . . . . . . 5.2 Deﬁnition of the DFT and its properties . . . . . . . . 5.2.1 Properties of the DFT . . . . . . . . . . . . . . 5.2.2 z-transforms ∗ . . . . . . . . . . . . . . . . . . . 5.3 Frequency description of wavelet analysis . . . . . . . 5.3.1 Low-pass and high-pass ﬁltering ∗ . . . . . . . . 5.4 Correlation and feature detection . . . . . . . . . . . . 5.4.1 DFT method of computing correlations . . . . 5.4.2 Proof of DFT eﬀect on correlation ∗ . . . . . . 5.4.3 Normalized correlations and feature detection ∗ 5.5 Object detection in 2D images . . . . . . . . . . . . . 5.6 Creating scaling signals and wavelets ∗ . . . . . . . . . 5.7 Gabor transforms and spectrograms . . . . . . . . . . 5.8 Musical analysis . . . . . . . . . . . . . . . . . . . . . 5.8.1 Analysis of Stravinsky’s Firebird Suite . . . . 5.8.2 Analysis of a Chinese folk song . . . . . . . . . 5.9 Inverting Gabor transforms . . . . . . . . . . . . . . . 5.10 Gabor transforms and denoising . . . . . . . . . . . . . 5.11 Notes and references . . . . . . . . . . . . . . . . . . . 5.12 Examples and exercises . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . .

167 168 169 170 171 173 174 178 180 181 183 183 185 188 192 195 197 199 201 203 206 210

6 Beyond wavelets 6.1 Wavelet packet transforms . . . . . . . . . . . . . 6.2 Wavelet packet transforms—applications . . . . . 6.2.1 Fingerprint compression . . . . . . . . . . 6.3 Continuous wavelet transforms . . . . . . . . . . 6.4 Gabor wavelets and speech analysis . . . . . . . . 6.4.1 Musical analysis: formants in song lyrics . 6.5 Percussion scalograms and musical rhythm . . . 6.5.1 Analysis of a complex percussive rhythm 6.5.2 Multiresolution Principle for rhythm . . . 6.6 Notes and references . . . . . . . . . . . . . . . . 6.6.1 Additional references . . . . . . . . . . . . 6.7 Examples and exercises . . . . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

223 223 225 228 228 232 236 237 241 241 241 242 246

A Projects A.1 Music . . . . . . . . . . . A.2 Noise removal from audio A.3 Wavelet image processing A.4 References . . . . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

255 255 256 256 257

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

B Selected exercise B.1 Introduction . B.2 Chapter 2 . . B.3 Chapter 3 . . B.4 Chapter 4 . . B.5 Chapter 5 . . B.6 Chapter 6 . .

solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

259 259 259 262 268 273 279

C Wavelet software 283 C.1 Installing the book’s software . . . . . . . . . . . . . . . . . . . 284 C.2 Other software . . . . . . . . . . . . . . . . . . . . . . . . . . . 284 C.3 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 285 Bibliography

287

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

Chapter 1

Overview

Thought is impossible without an image. —Aristotle

We provide a brief, non-mathematical, summary of what wavelet analysis is all about. Examples are discussed from imaging and audio. In the Notes and References section we indicate where some good material on the history of wavelet analysis can be found.

1.1

What is wavelet analysis?

Wavelet analysis decomposes sounds and images into component waves of varying durations, called wavelets. These wavelets are localized vibrations of a sound signal, or localized variations of detail in an image. They can be used for a wide variety of fundamental signal processing tasks, such as compression, removing noise, or enhancing a recorded sound or image in various ways. In this section we give a couple of illustrations of what wavelet analysis involves, chosen from the ﬁelds of imaging and audio. In line with the dictum of Aristotle quoted above, we shall ﬁrst consider the ﬁeld of imaging. The heart of wavelet analysis is multiresolution analysis. Multiresolution analysis is the decomposition of a signal (such as an image) into subsignals (subimages) of diﬀerent size resolution levels. As an example, consider how the artist Ray Freer describes his process of painting: My approach to watercolor painting is very similar to that with oil or pastel. I start with a broad description of the essential masses and then © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

1. Overview

2

progress through several stages of reﬁnement until I have achieved the ﬁnal result.1

To see what Freer means, look at Figure 1.1. In Figure 1.1(a) we show an initial image, corresponding to the “broad description of the essential masses.” Figure 1.1(d) shows the details that need to be added to image (a) in order to obtain the ﬁrst reﬁnement shown in Figure 1.1(b). These details consist of regions of black where image (a) is to be darkened, and white where image (a) is to be brightened, and a dull gray background where image (a) is left unchanged. If we view these regions as a surface of water (the white parts as higher, the black parts as lower, above the undisturbed gray surface) then the white and black regions in (d) are localized waves, wavelets, that represent the changes to be made to image (a). Similarly, in image (e) we show the wavelets representing the changes to be made to image (b) to produce the ﬁnal image (c). Notice that the wavelets in (e) are ﬁner in detail, smaller scale, than the wavelets in (d). In wavelet analysis, we begin with the reverse of the process shown in Figure 1.1. We start with the ﬁnal image and determine the ﬁnest scale wavelets contained in it, and remove those to give a less detailed image. We then determine the next ﬁnest scale wavelets contained in that image and remove those, and so on, until we arrive at an image consisting of “broad masses,” a very low resolution image. This process is called the wavelet transform of the image. The wavelet transform thus provides us with the recipe for synthesizing the image as we described above. What advantages do we gain from this wavelet transform process? To see what we gain in terms of the painting example just discussed, we quote from Enrico Coen’s excellent book The Art of Genes: Why doesn’t the artist go directly to the detailed last stage, circumventing the other steps? What purpose do the previous stages serve? One problem with going too quickly into painting details is that they can be easily misplaced. If you start by painting the eyes in great detail, you might ﬁnd later on that they are not quite in the right place in relation to the other features. You would have to go back and start all over again. By ﬁrst giving the overall layout and gradually narrowing down to the details, the correct relations between features can be established gradually at each scale, from the overall location of the face, to the position of the eyes, to the details within each eye. To help them achieve this, some artists screw up their eyes when looking at the subject during the early stages of a painting, blurring their vision on purpose so as to more easily see the broad arrangements and avoid being distracted by details.2

The last sentence quoted gives us one clue as to how wavelet analysis removes random noise (random clutter) from images. If you “screw up” your eyes to 1 From 2 From

[1], page 132. [1], page 133.

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

1.1 What is wavelet analysis?

(a)

(b)

(d)

3

(c)

(e)

FIGURE 1.1 Drawing a portrait. Figure (a) is the first, rough draft. Figure (d) shows details that are added to (a) to create the second draft (b). In (d), points are blacker if the corresponding points in (a) are to be darkened. They are whiter if the corresponding points in (a) are to be brightened. And they are a dull gray if no change is needed in (a). Similarly, in (e) we show details that are needed to change (b) into the third, and final, draft (c).

blur your vision when looking at one of the noisy images in Chapter 4 [say Figure 4.17(b) on page 137, or Figure 4.18(a) on page 138], you will see that the noise is removed from the blurred image. Of course, there is much more to wavelet-based denoising than that. We will show in Chapter 4 how to add back in details to obtain sharply focused images with the noise removed. But what does this description of wavelet analysis have to do with compression of images? We can give some idea in the following way. Notice that the wavelets in Figure 1.1(d) are of a signiﬁcantly larger scale than the ones in Figure 1.1(e). This corresponds to relatively fewer brushstrokes (with a broader brush) needed to reﬁne the image in (a) to produce the image in (b). Likewise, to produce the ﬁrst draft image (a), even fewer and broader brushstrokes are needed. Because we need progressively fewer brushstrokes at these lower resolutions, the image can be compressed by specifying how to produce the lower resolution, ﬁrst draft, image and then how to successively reﬁne it (rather than trying to specify how to paint in all the ﬁne details in the original image). We will see in Chapter 4 that this, in essence, is what wavelet transform compression is doing. While 2D images are easier to see, we will begin our mathematical presentation in the next chapter with 1D signals because the mathematics is simpler © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

1. Overview

4

in 1D. The overall process, however, is similar to what we described above for images. The wavelet transform speciﬁes how to remove the ﬁnest details from a sound signal (the shortest duration wavelets) to produce a lower resolution signal. But what does “lower resolution” mean for a sound signal? It means a signal where the oscillations are of longer duration in comparison to the short duration wavelets subtracted out from the original. Steven Pinker describes this very well in terms of musical sounds: It [musical analysis] dissects the melody into essential parts and ornaments. The ornaments are stripped oﬀ and the essential parts further dissected into even more essential parts and ornaments on them. . . we sense it when we recognize variations of a piece in classical music or jazz. The skeleton of the melody is conserved while the ornaments diﬀer from variation to variation.3

The ﬁrst two sentences of this quotation are an excellent, non-mathematical, summary of what the wavelet transform is doing to 1D signals. Keep it in mind in the next two chapters when we discuss 1D wavelet transforms. We will begin with the Haar wavelet transform in Chapter 2. This transform provides the simplest way to decompose a 1D signal by stripping oﬀ details (called fluctuations) to obtain a skeleton (called a trend ). After setting the stage with Haar wavelets, we then turn to the Daubechies wavelet transforms in Chapter 3. These transforms provide a powerful approach to signal compression and signal denoising, and a deep understanding of signal synthesis.

1.2

Notes and references

For those readers interested in the history of wavelet analysis, a good place to start is the article by Burke [3], which has been expanded into a book [4]. There is also some history, of a more sophisticated mathematical nature, in the books by Meyer [5] and Daubechies [6]. 1. E. Coen. (1999). The Art of Genes. Oxford University Press, Oxford, UK. 2. S. Pinker. (1997). How the Mind Works. Norton, NY. 3. B. Burke. (1994). The mathematical microscope: waves, wavelets, and beyond. In A Positron Named Priscilla, Scientific Discovery at the Frontier, 196–235, Nat. Academy Press. Edited by M. Bartusiak. 4. B. Burke Hubbard. (1998). The World According to Wavelets, Second Edition. AK Peters, Wellesley, MA. 5. Y. Meyer. (1993). Wavelets. Algorithms and Applications. SIAM, Philadelphia, PA. 6. I. Daubechies. (1992). Ten Lectures on Wavelets. SIAM, Philadelphia, PA.

3 From

[2], page 533.

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

Chapter 2

Haar wavelets

The purpose of computing is insight, not numbers. —Richard W. Hamming

The purpose of computing is insight, not pictures. —Lloyd N. Trefethen1

A Haar wavelet is the simplest type of wavelet. In discrete form, Haar wavelets are related to a mathematical operation called the Haar transform. The Haar transform serves as a prototype for all other wavelet transforms. Studying the Haar transform in detail will provide a good foundation for understanding the more sophisticated wavelet transforms which we shall describe in the next chapter. In this chapter we shall describe how the Haar transform can be used for compressing audio signals and for removing noise. Our discussion of these applications will set the stage for the more powerful wavelet transforms to come and their applications to these same problems. One distinctive feature that the Haar transform enjoys is that it lends itself easily to simple hand calculations. We shall illustrate many concepts by both simple hand calculations and more involved computer computations.

1 Hamming’s

quote is from [1]. Trefethen’s quote is from [2].

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

6

2.1

The Haar transform

In this section we shall introduce the basic notions connected with the Haar transform, which we shall examine in more detail in later sections. First, we need to deﬁne the type of signals that we shall be analyzing. We shall be working extensively with discrete signals. A discrete signal is a function of time with values occurring at discrete instants. Generally we shall express a discrete signal in the form f = (f1 , f2 , . . . , fN ), where N is a positive even integer2 which we shall refer to as the length of f . The values of f are the N real numbers f1 , f2 , . . . , fN . These values are typically measured values of an analog signal g, measured at the time values t = t1 , t2 , . . . , tN . That is, the values of f are f1 = g(t1 ), f2 = g(t2 ), . . . , fN = g (tN ) .

(2.1)

For simplicity, we shall assume that the increment of time that separates each pair of successive time values is always the same. We shall use the phrase equally spaced sample values, or just sample values, when the discrete signal has its values deﬁned in this way. An important example of sample values is the set of data values stored in a computer audio ﬁle, such as a .wav ﬁle. Another example is the sound intensity values recorded on a compact disc. A non-audio example is a digitized electrocardiogram. Like all wavelet transforms, the Haar transform decomposes a discrete signal into two subsignals of half its length. One subsignal is a running average or trend; the other subsignal is a running diﬀerence or ﬂuctuation. Let’s begin with the trend. The ﬁrst trend, a1 = (a1 , a2 , . . . , aN/2 ), for the signal f is computed by taking a running average in the following way. Its ﬁrst average of the ﬁrst pair of value, a1 , is computed by taking the √ √ values of f : (f1 +f2 )/2; and then multiplying it by 2. Thus, a1 = (f1 +f2 )/ 2. Similarly, of the next pair of values its next value a2 is computed by taking the average √ √ of f : (f3 + f4 )/2; and then multiplying it by 2. Hence, a2 = (f3 + f4 )/ 2. Continuing in this way, all of the values of a1 are produced by taking averages √ of successive pairs of values of f , and then multiplying these averages by 2. A formula for the values of a1 is am =

f2m−1 + f2m √ , 2

(2.2)

for m = 1, 2, 3, . . . , N/2. For example, suppose f is deﬁned by eight values, say f = (4, 6, 10, 12, 8, 6, 5, 5). 2 If

N is odd, then a zero can be appended to the signal to make it even.

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.1 The Haar transform

7

The average of its ﬁrst two values is 5, the average of the next two values is 11, the average of the next two values is√7, and the average of the last two values is 5. Multiplying these averages by 2 , we obtain the ﬁrst trend subsignal √ √ √ √ a1 = (5 2, 11 2, 7 2, 5 2). √ You might ask: Why perform the extra step of multiplying by 2 ? Why not just take averages? These questions √ will be answered in the next section, when we show that multiplication by 2 is needed in order to ensure that the Haar transform preserves the energy of a signal. The other subsignal is called the ﬁrst ﬂuctuation. The ﬁrst ﬂuctuation of the signal f , which is denoted by d1 = (d1 , d2 , . . . , dN/2 ), is computed by taking a running diﬀerence in the following way. Its ﬁrst value, d1 , is found by taking half the√diﬀerence of the ﬁrst pair√of values of f : (f1 − f2 )/2; and multiplying it by 2. That is, d1 = (f1 −f2 )/ 2. Likewise, its next value d2 is found by taking half the √ f : (f3 − f4 )/2; √diﬀerence of the next pair of values of and multiplying it by 2. In other words, d2 = (f3 − f4 )/ 2. Continuing in this way, all of the values of d1 are produced according to the following formula: f2m−1 − f2m √ , (2.3) dm = 2 for m = 1, 2, 3, . . . , N/2. For example, for the signal considered above, f = (4, 6, 10, 12, 8, 6, 5, 5), its ﬁrst ﬂuctuation is

2.1.1

√ √ √ d1 = (− 2, − 2, 2, 0).

Haar transform, 1-level

The Haar transform is performed in several stages, or levels. The ﬁrst level is the mapping H1 deﬁned by 1 (a1 | d1 ) f −→

H

(2.4)

from a discrete signal f to its ﬁrst trend a1 and ﬁrst ﬂuctuation d1 . For example, we showed above that √ √ √ √ √ √ √ H1 (4, 6, 10, 12, 8, 6, 5, 5) −→ (5 2, 11 2, 7 2, 5 2 | − 2, − 2, 2, 0).

(2.5)

The mapping H1 in (2.4) has an inverse. Its inverse maps the transform signal (a1 | d1 ) back to the signal f , via the following formula: aN/2 + dN/2 aN/2 − dN/2 a1 + d1 a1 − d1 √ √ √ , √ ,..., , . (2.6) f= 2 2 2 2 © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

8

√ √ √ In other words,√f1 = (a1 + d1 )/ 2, f2 = (a1 − d1 )/ 2, f3 = (a2 + d2 )/ 2, f4 = (a2 − d2 )/ 2, and so on. The reader may wish to check that the formulas in (2.6) when applied to √ √ √ √ √ √ √ (5 2, 11 2, 7 2, 5 2 | − 2, − 2, 2, 0) produce (4, 6, 10, 12, 8, 6, 5, 5) thereby inverting the transform in (2.5). Let’s now consider what advantages accrue from performing the Haar transformation. These advantages will be described in more detail later in this chapter, but some basic notions can be introduced now. All of these advantages stem from the following cardinal feature of the Haar transform (a feature that will be even more prominent for the Daubechies transforms described in the next chapter): Small Fluctuations Feature. The magnitudes of the values of the ﬂuctuation subsignal are often signiﬁcantly smaller than the magnitudes of the values of the original signal. For instance, for the signal f = (4, 6, 10, 12, 8, 6, 5, 5) considered above, its eight values have an √ average of 7. On the other hand, for its ﬁrst √ magnitude √ 1 ﬂuctuation d = (− 2, − 2, 2, 0), the average of its four magnitudes is √ 0.75 2. In this case, the magnitudes of the ﬂuctuation’s values are an average of 6.6 times smaller than the magnitudes of the original signal’s values. For a second example, consider the signal shown in Figure 2.1(a). This signal was generated from 1024 sample values of the function g(x) = 20x2 (1 − x)4 cos 12πx over the interval [0, 1). In Figure 2.1(b) we show a graph of the 1-level Haar transform of this signal. The trend subsignal is graphed on the left half, over the interval [0, 0.5), and the ﬂuctuation subsignal is graphed on the right half, over the interval [0.5, 1). It is clear that a large percentage of the ﬂuctuation’s values are close to 0 in magnitude, another instance of the Small Fluctuations Feature. Notice also that the trend subsignal looks like the original signal, √ although shrunk by half in length and expanded by a factor of 2 vertically. The reason that the Small Fluctuations Feature is generally true is that typically we are dealing with signals sampled from a continuous analog signal g with a very short time increment between the samples. In other words, the equations in (2.1) hold with a small value of the time increment h = tk+1 − tk for each k = 1, 2, . . . , N − 1. If the time increment is small enough, then successive values f2m−1 = g (t2m−1 ) and f2m = g (t2m ) of the signal f will be close to each other due to the continuity of g. Consequently, the ﬂuctuation values for the Haar transform satisfy dm =

g (t2m−1 ) − g (t2m ) √ ≈ 0. 2

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.2 Conservation and compaction of energy

0

0.25

0.5

0.75

9

1

1

0.5

0.5

0

0

−0.5

−0.5

−1 1

(a)

0

0.25

0.5

0.75

−1 1

(b)

FIGURE 2.1 (a) Signal, (b) Haar transform, 1-level.

This explains why the Small Fluctuations Feature is generally true for the Haar transform. A similar analysis shows why the trend subsignal has a graph that is similar in appearance to the ﬁrst trend. If g is continuous and the time increment is very small, then g(t2m−1 ) and g(t2m ) will be close to each other. Expressing this fact as an approximation, g(t2m−1 ) ≈ g(t2m ), we obtain the following approximation for each value am of the trend subsignal √ am ≈ 2 g(t2m ). This equation shows that a1 is approximately the same as sample values of √ 2 g(x) for x = t2 , t4 , . . . , tN . In other words, it shows that the graph of the ﬁrst trend subsignal is similar in appearance to the graph of g, as we pointed out above in regard to the signal in Figure 2.1(a). We shall examine these points in more detail in the next chapter. One of the reasons that the Small Fluctuations Feature is important is that it has applications to signal compression. By compressing a signal we mean transmitting its values, or approximations of its values, by using a smaller number of bits. For example, if we were only to transmit the trend subsignal for the signal shown in Figure 2.1(a) and then perform Haar transform inversion (treating the ﬂuctuation’s values as all zeros), then we would obtain an approximation of the original signal. Since the length of the trend subsignal is half the length of the original signal, this would achieve 50% compression. We shall discuss compression in more detail in Section 2.5. Once we have performed a 1-level Haar transform, then it is easy to perform multiple-level Haar transforms. We shall discuss this in the next section.

2.2

Conservation and compaction of energy

In the previous section we deﬁned the 1-level Haar transform. In this section we shall discuss its two most important properties: (1) It conserves the ener© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

10

gies of signals; (2) It performs a compaction of the energy of signals. We shall also complete our deﬁnition of the Haar transform by showing how to extend its deﬁnition to multiple levels.

2.2.1

Conservation of energy

An important property of the Haar transform is that it conserves the energies of signals. By the energy of a signal f we mean the sum of the squares of its values. That is, the energy Ef of a signal f is deﬁned by 2 Ef = f12 + f22 + · · · + fN .

(2.7)

We shall provide some explanation for why we use the term energy for this quantity Ef in a moment. First, however, let’s look at an example of calculating energy. Suppose f = (4, 6, 10, 12, 8, 6, 5, 5) is the signal considered in Section 2.1. Then Ef is calculated as follows: Ef = 42 + 62 + · · · + 52 = 446. Furthermore, using the values for its 1-level Haar transform √ √ √ √ √ √ √ (a1 | d1 ) = (5 2, 11 2, 7 2, 5 2 | − 2, − 2, 2, 0), we ﬁnd that E(a1 | d1 ) = 25 · 2 + 121 · 2 + · · · + 2 + 0 = 446 as well. Thus the 1-level Haar transform has kept the energy constant. In fact, this is true in general: Conservation of energy. The 1-level Haar transform conserves energy. In other words, E(a1 | d1 ) = Ef for every signal f . We will explain why this conservation of energy property is true for all signals at the end of this section. Before we go any further, we should say something about why we use the term energy for the quantity Ef . The reason is that sums of squares frequently appear in physics when various types of energy are calculated. For instance, if a particle of mass m has a velocity of v = (v1 , v2 , v3 ), then its kinetic energy is (m/2)(v12 + v22 + v32 ). Hence its kinetic energy is proportional to v12 + v22 + v32 = Ev . Ignoring the constant of proportionality, m/2, we obtain the quantity Ev which we call the energy of v. While conservation of energy is certainly an important property, it is even more important to consider how the Haar transform redistributes the energy in a signal by compressing most of the energy into the trend subsignal. For example, for the signal √ f = (4, √ 6, 10, √ 12,√8, 6, 5, 5) we found in Section 2.1 that its trend a1 equals (5 2, 11 2, 7 2, 5 2). Therefore, the energy of a1 is Ea1 = 25 · 2 + 121 · 2 + 49 · 2 + 25 · 2 = 440. © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.2 Conservation and compaction of energy

11

√ √ √ On the other hand, the ﬂuctuation d1 is (− 2, − 2, 2, 0), which has energy Ed1 = 2 + 2 + 2 + 0 = 6. Thus the energy of the trend a1 accounts for 440/446 = 98.7% of the total energy of the signal. In other words, the 1-level Haar transform has redistributed the energy of f so that over 98% is concentrated into the subsignal a1 which is half the length of f . For obvious reasons, this is called compaction of energy. As another example, consider the signal f graphed in Figure 2.1(a) and its 1-level Haar transform shown in Figure 2.1(b). In this case, we ﬁnd that the energy of the signal f is 127.308 while the energy of its ﬁrst trend a1 is 127.305. Thus 99.998% of the total energy is compacted into the half-length subsignal a1 . By examining the graph in Figure 2.1(b) it is easy to see why such a phenomenal energy compaction has occurred; the values of the ﬂuctuation d1 are so small, relative to the much larger values of the trend a1 , that its energy Ed1 contributes only a small fraction of the total energy Ea1 + Ed1 . These two examples illustrate the following general principle: Compaction of energy. The energy of the trend subsignal a1 accounts for a large percentage of the energy of the transformed signal (a1 | d1 ). Compaction of energy will occur whenever the magnitudes of the ﬂuctuation’s values are signiﬁcantly smaller than the trend’s values (recall the Small Fluctuations Feature from the last section). In Section 2.5, we shall describe how compaction of energy provides a framework for applying the Haar transform to compress signals. We now turn to a discussion of how the Haar transform can be extended to multiple levels, thereby increasing the energy compaction of signals.

2.2.2

Haar transform, multiple levels

Once we have performed a 1-level Haar transform, then it is easy to repeat the process and perform multiple level Haar transforms. After performing a 1-level Haar transform of a signal f we obtain a ﬁrst trend a1 and a ﬁrst ﬂuctuation d1 . The second level of a Haar transform is then performed by computing a second trend a2 and a second ﬂuctuation d2 for the ﬁrst trend a1 only. For example, if f = (4, 6, 10, 12, 8, 6, 5, 5)√is the above, √ signal √ considered √ then we found that its ﬁrst trend is a1 = (5 2, 11 2, 7 2, 5 2). To get the 1 second trend a2 we apply Formula (2.2)√to the√values √ of√ a . That is, we 1 add successive pairs of values of a = (5 2, 11 2, 7 2, 5 2) and divide by √ 2 obtaining a2 = (16, 12). To get √ the second ﬂuctuation d2 , we subtract √ √ √ √ 1 successive pairs of values of a = (5 2, 11 2, 7 2, 5 2) and divide by 2 obtaining a2 = (−6, 2). Thus the 2-level Haar transform of f is the signal √ √ √ (a2 | d2 | d1 ) = (16, 12 | −6, 2 | − 2, − 2, 2, 0). © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

12

For this signal f , a 3-level Haar transform can also be done, and the result is √ √ √ √ √ (a3 | d3 | d2 | d1 ) = (14 2 | 2 2 | −6, 2 | − 2, − 2, 2, 0). It is interesting to calculate the energy compaction that has occurred with the 2-level and 3-level Haar transforms that we just computed. First, we know that E(a2 | d2 | d1 ) = 446 because of conservation of energy. Second, we compute that Ea2 = 400. Thus the 2-level Haar transformed signal (a2 | d2 | d1 ) has almost 90% of the total energy of f contained in the second trend a2 which is 1/4 of the length of f . This is a further compaction, or localization, of the energy of f . Furthermore, Ea3 = 392; thus a3 contains 87.89% of the total energy of f . This is even further compaction; the 3-level Haar transform (a3 | d3 | d2 | d1 ) has almost 88% of the total energy of f contained in the third trend a3 which is 1/8 the length of f . For those readers who are familiar with Quantum Theory, there is an interesting phenomenon here that is worth noting. By Heisenberg’s Uncertainty Principle, it is impossible to localize a ﬁxed amount of energy into an arbitrarily small time interval. This provides an explanation for why the energy percentage dropped from 98% to 90% when the second-level Haar transform was computed, and from 90% to 88% when the third-level Haar transform was computed. When we attempt to squeeze the energy into ever smaller time intervals, it is inevitable that some energy leaks out. As another example of how the Haar transform redistributes and localizes the energy in a signal, consider the graphs shown in Figure 2.2. In Figure 2.2(a) we show a signal, and in Figure 2.2(b) we show the 2-level Haar transform of this signal. In Figures 2.2(c) and (d) we show the respective cumulative energy proﬁles of these two signals. By the cumulative energy proﬁle of a signal f we mean the signal deﬁned by

f12 f12 + f22 f12 + f22 + f32 , , , ...,1 . Ef Ef Ef

The cumulative energy proﬁle of f thus provides a summary of the accumulation of energy in the signal as time proceeds. As can be seen from comparing the two proﬁles in Figure 2.2, the 2-level Haar transform has redistributed and localized the energy of the original signal.

2.2.3

Justification of conservation of energy

We close this section with a brief justiﬁcation of the conservation of energy property of the Haar transform. First, we observe that the terms a21 and d21 © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.2 Conservation and compaction of energy

0

0.5

1

1.5

1.5

1.5

0.75

0.75

0

0

−0.75

−0.75

−1.5 2

(a)

0

13

0

0.5

1

1.5

−1.5 2

(b)

0.5

1

1.5

(c)

1.5

1.5

1

1

0.5

0.5

0

0

−0.5 2

0

0.5

1

1.5

−0.5 2

(d)

FIGURE 2.2 (a) Signal. (b) 2-level Haar transform of signal. (c) Cumulative energy profile of signal. (d) Cumulative energy profile of 2-level Haar transform.

in the formula E(a1 | d1 ) = a21 + · · · + a2N/2 + d21 + · · · + d2N/2 add up as follows: 2 2 f1 + f2 f1 − f2 √ √ + 2 2 2 2 f 2 − 2f1 f2 + f22 f + 2f1 f2 + f2 + 1 = 1 2 2 = f12 + f22 .

a21 + d21 =

2 2 + f2m for all other values of m. Therefore, by Similarly, a2m + d2m = f2m−1 adding a2m and d2m successively for each m, we ﬁnd that 2 . a21 + · · · + a2N/2 + d21 + · · · + d2N/2 = f12 + · · · + fN

In other words, E(a1 | d1 ) = Ef , which justiﬁes the conservation of energy property. © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

14

2.3

Haar wavelets

In this section we discuss the simplest wavelets, the Haar wavelets. This material will set the stage for the more sophisticated Daubechies wavelets described in the next chapter. We begin by discussing the 1-level Haar wavelets. These wavelets are deﬁned as 1 −1 1 W1 = √ , √ , 0, 0, . . . , 0 2 2 1 −1 1 W2 = 0, 0, √ , √ , 0, 0, . . . , 0 2 2 −1 1 1 √ √ W3 = 0, 0, 0, 0 , , 0, 0, . . . , 0 2 2 .. . 1 −1 1 WN/2 = 0, 0, . . . , 0, √ , √ . (2.8) 2 2 These 1-level Haar wavelets have a number of interesting properties. First, they each have an energy of 1. Second, √ they each consist of a rapid ﬂuctuation between just two non-zero values, ±1/ 2, with an average value of 0. Hence the name wavelet. Finally, they all are very similar to each other in that they are each a translation in time by an even number of time-units of the ﬁrst Haar wavelet W11 . The second Haar wavelet W21 is a translation forward in time by two units of W11 , and W31 is a translation forward in time by four units of W11 , and so on. The reason for introducing the 1-level Haar wavelets is that we can express the 1-level ﬂuctuation subsignal in a simpler form by making use of scalar products with these wavelets. The scalar product is a fundamental operation on two signals, and is deﬁned as follows. Scalar product: The scalar product f · g of the signals f = (f1 , f2 , . . . , fN ) and g = (g1 , g2 , . . . , gN ) is deﬁned by f · g = f1 g1 + f2 g2 + · · · + fN gN .

(2.9)

Using the 1-level Haar wavelets, we can express the values for the ﬁrst ﬂuctuation subsignal d1 as scalar products. For example, d1 =

f1 − f2 √ = f · W11 . 2

Similarly, d2 = f · W21 , and so on. We can summarize Formula (2.3) in terms of scalar products with the 1-level Haar wavelets: 1 dm = f · Wm

(2.10)

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.3 Haar wavelets

15

for m = 1, 2, . . . , N/2. We can also use the idea of scalar products to restate the Small Fluctuations Feature from Section 2.1 in a more precise form. If we say that the support of each Haar wavelet is the set of two time-indices where the wavelet is non-zero, then we have the following more precise version of the Small Fluctuations Feature: Property 1. If a signal f is (approximately) constant over the support of a 1-level Haar wavelet Wk1 , then the ﬂuctuation value dk = f · Wk1 is (approximately) zero. This property will be considerably strengthened in the next chapter. Note: From now on, we shall refer to the set of time-indices m where fm = 0 as the support of a signal f . For example, the support of the signal f = (0, 0, 5, 7, −2, 0, 2, 0) consists of the indices 3, 4, 5, and 7; while the support of the Haar wavelet W21 consists of the indices 3, 4. We can also express the 1-level trend values as scalar products with certain elementary signals. These elementary signals are called 1-level Haar scaling signals, and they are deﬁned as 1 1 = √ , √ , 0, 0, . . . , 0 2 2 1 1 1 V2 = 0, 0, √ , √ , 0, 0, . . . , 0 2 2 .. . 1 1 1 = 0, 0, . . . , 0, √ , √ . VN/2 2 2 V11

(2.11)

Using these Haar scaling signals, the values a1 , . . . , aN/2 for the ﬁrst trend are expressed as scalar products: 1 am = f · Vm

(2.12)

for m = 1, 2, . . . , N/2. The Haar scaling signals are quite similar to the Haar wavelets. They all have energy 1 and have a support consisting of just two consecutive timeindices. In fact, they are all translates by an even multiple of time-units of the ﬁrst scaling signal V11 . Unlike the Haar wavelets, however, the average values of the Haar √scaling signals are not zero. In fact, they each have an average value of 1/ 2. The ideas discussed above extend to every level. For simplicity, we restrict our discussion to the second level. The 2-level Haar scaling signals are deﬁned © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

16 by

1 1 1 1 , , , , 0, 0 . . . , 0 2 2 2 2 1 1 1 1 V22 = 0, 0, 0, 0, , , , , 0, 0, . . . , 0 2 2 2 2 .. . 1 1 1 1 2 . VN/4 = 0, 0, . . . , 0, , , , 2 2 2 2 V12 =

(2.13)

These scaling signals are all translations by multiples of four time-units of the ﬁrst scaling signal V12 , and they all have energy 1 and average value 1/2. Furthermore, the values of the 2-level trend a2 are scalar products of these scaling signals with the signal f . That is, a2 satisﬁes 2 . (2.14) a2 = f · V12 , f · V22 , . . . , f · VN/4 Likewise, the 2-level Haar wavelets are deﬁned by 1 1 −1 −1 W12 = , , , , 0, 0 . . . , 0 2 2 2 2 1 1 −1 −1 W22 = 0, 0, 0, 0, , , , , 0, 0, . . . , 0 2 2 2 2 .. . 1 1 −1 −1 2 , . = 0, 0, . . . , 0, , , WN/4 2 2 2 2

(2.15)

These wavelets all have supports of length 4, since they are all translations by multiples of four time-units of the ﬁrst wavelet W12 . They also all have energy 1 and average value 0. Using scalar products, the 2-level ﬂuctuation d2 satisﬁes 2 . (2.16) d2 = f · W12 , f · W22 , . . . , f · WN/4

2.4

Multiresolution analysis

In the previous section we discussed how the Haar transform can be described using scalar products with scaling signals and wavelets. In this section we discuss how the inverse Haar transform can also be described in terms of these same elementary signals. This discussion will show how discrete signals are synthesized by beginning with a very low resolution signal and successively adding on details to create higher resolution versions, ending with a complete synthesis of the signal at the ﬁnest resolution. This is known as multiresolution analysis (MRA). MRA is the heart of wavelet analysis. © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.4 Multiresolution analysis

17

In order to make these ideas precise, we must ﬁrst discuss some elementary operations that can be performed on signals. Given two signals f = (f1 , f2 , . . . , fN ) and g = (g1 , g2 , . . . , gN ), we can perform the following algebraic operations: Addition and Subtraction: The sum f + g of the signals f and g is deﬁned by adding their values: f + g = (f1 + g1 , f2 + g2 , . . . , fN + gN ).

(2.17)

Their diﬀerence f − g is deﬁned by subtracting their values: f − g = (f1 − g1 , f2 − g2 , . . . , fN − gN ).

(2.18)

Constant multiple: A signal f is multiplied by a constant c by multiplying each of its values by c. That is, c f = (cf1 , cf2 , . . . , cfN ).

(2.19)

For example, by repeatedly applying the addition operation, we can express a signal f = (f1 , f2 , . . . , fN ) as follows: f = (f1 , 0, 0, . . . , 0) + (0, f2 , 0, 0, . . . , 0) + · · · + (0, 0, . . . , 0, fN ). Then, by applying the constant multiple operation to each of the signals on the right side of this last equation, we obtain f = f1 (1, 0, 0, . . . , 0) + f2 (0, 1, 0, 0, . . . , 0) + · · · + fN (0, 0, . . . , 0, 1). This formula is a very natural one; it amounts to expressing f as a sum of its individual values at each discrete instant of time. 0 as If we deﬁne the elementary signals V10 , V20 , . . . , VN V10 = (1, 0, 0, . . . , 0) V20 = (0, 1, 0, 0, . . . , 0) .. . 0 VN = (0, 0, . . . , 0, 1)

(2.20)

then the last formula for f can be rewritten as 0 . f = f1 V10 + f2 V20 + · · · + fN VN

(2.21)

Formula (2.21) is called the natural expansion of a signal f in terms of the 0 . We shall now show that the Haar natural basis of signals V10 , V20 , . . . , VN MRA involves expressing f as a sum of constant multiples of a diﬀerent basis set of elementary signals: the Haar wavelets and scaling signals deﬁned in the previous section. © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

18

In the previous section, we showed how to express the 1-level Haar transform in terms of wavelets and scaling signals. It is also possible to express the inverse of the 1-level Haar transform in terms of these same elementary signals. This leads to the ﬁrst level of the Haar MRA. To deﬁne this ﬁrst level Haar MRA we make use of (2.6) to express a signal f as aN/2 aN/2 a1 a1 a2 a2 f = √ , √ , √ , √ ,..., √ , √ 2 2 2 2 2 2 dN/2 −dN/2 d1 −d1 d2 −d2 . + √ , √ , √ , √ ,..., √ , √ 2 2 2 2 2 2 This formula shows that the signal f can be expressed as the sum of two signals that we shall call the ﬁrst averaged signal and the ﬁrst detail signal. That is, we have f = A1 + D1 (2.22) where the signal A1 is called the ﬁrst averaged signal and is deﬁned by aN/2 aN/2 a1 a1 a2 a2 1 A = √ , √ , √ , √ ,..., √ , √ (2.23) 2 2 2 2 2 2 and the signal D1 is called the ﬁrst detail signal and is deﬁned by dN/2 −dN/2 d1 −d1 d2 −d2 1 D = √ , √ , √ , √ ,..., √ , √ . 2 2 2 2 2 2

(2.24)

Using Haar scaling signals and wavelets, and using the basic algebraic operations with signals, the averaged and detail signals are expressible as 1 , A1 = a1 V11 + a2 V21 + · · · + aN/2 VN/2

(2.25a)

1 D1 = d1 W11 + d2 W21 + · · · + dN/2 WN/2 .

(2.25b)

Applying the scalar product formulas for the coeﬃcients in Equations (2.10) and (2.12), we can rewrite these last two formulas as follows 1 1 )VN/2 , A1 = (f · V11 )V11 + (f · V21 )V21 + · · · + (f · VN/2 1 1 )WN/2 . D1 = (f · W11 )W11 + (f · W21 )W21 + · · · + (f · WN/2

These formulas show that the averaged signal is a combination of Haar scaling signals, with the values of the ﬁrst trend subsignal as coeﬃcients; and that the detail signal is a combination of Haar wavelets, with the values of the ﬁrst ﬂuctuation subsignal as coeﬃcients. As an example of these ideas, consider the signal f = (4, 6, 10, 12, 8, 6, 5, 5). © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.4 Multiresolution analysis

19

In Section 2.1 we found that its ﬁrst trend subsignal was √ √ √ √ a1 = (5 2, 11 2, 7 2, 5 2). Applying Formula (2.23), the averaged signal is A1 = (5, 5, 11, 11, 7, 7, 5, 5).

(2.26)

Notice how the ﬁrst averaged signal consists of the repeated average values 5, 5, and 11, 11, and 7, 7, and 5, 5 about which the values of f ﬂuctuate. Using Formula (2.25a), the ﬁrst averaged signal can also be expressed in terms of Haar scaling signals as √ √ √ √ A1 = 5 2 V11 + 11 2 V21 + 7 2 V31 + 5 2 V41 . Comparing these last two equations we can see that the positions of the repeated averages correspond precisely with the supports of the scaling signals. We also √ found √ √in Section 2.1 that the ﬁrst ﬂuctuation signal for f was d1 = (− 2, − 2, 2, 0). Formula (2.24) then yields D1 = (−1, 1, −1, 1, 1, −1, 0, 0). Thus, using the result for A1 computed above, we have f = (5, 5, 11, 11, 7, 7, 5, 5) + (−1, 1, −1, 1, 1, −1, 0, 0). This equation illustrates the basic idea of MRA. The signal f is expressed as a sum of a lower resolution, or averaged, signal (5, 5, 11, 11, 7, 7, 5, 5) added with a signal (−1, 1, −1, 1, 1, −1, 0, 0) made up of ﬂuctuations or details. These ﬂuctuations provide the added details necessary to produce the full resolution signal f . For this example, using Formula (2.25b), the ﬁrst detail signal can also be expressed in terms of Haar wavelets as √ √ √ D1 = − 2 W11 − 2 W21 + 2 W31 + 0 W41 . This formula shows that the values of D1 occur in successive pairs of rapidly ﬂuctuating values positioned at the supports of the Haar wavelets.

2.4.1

Multiresolution analysis, multiple levels

In the discussion above, we described the ﬁrst level of the Haar MRA of a signal. This idea can be extended to further levels, as many levels as the number of times that the signal length can be divided by 2. The second level of a MRA of a signal f involves expressing f as f = A 2 + D2 + D1 .

(2.27)

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

20

Here A2 is the second averaged signal and D2 is the second detail signal. Comparing Formulas (2.22) and (2.27) we see that A1 = A2 + D2 .

(2.28)

This formula expresses the fact that computing the second averaged signal A2 and second detail signal D2 simply consists of performing a ﬁrst level MRA of the signal A1 . Because of this, it follows that the second level averaged signal A2 satisﬁes 2 2 )VN/4 A2 = (f · V12 )V12 + (f · V22 )V22 + · · · + (f · VN/4

and the second level detail signal D2 satisﬁes 2 2 )WN/4 . D2 = (f · W12 )W12 + (f · W22 )W22 + · · · + (f · WN/4

For example, if f = (4, 6, 10, 12, 8, 6, 5, 5), then we found in Section 2.2 that a2 = (16, 12). Therefore 1 1 1 1 1 1 1 1 , , , , 0, 0, 0, 0 + 12 0, 0, 0, 0, , , , A2 = 16 2 2 2 2 2 2 2 2 = (8, 8, 8, 8, 6, 6, 6, 6).

(2.29)

It is interesting to compare the equations in (2.26) and (2.29). The second averaged signal A2 has values created from averages that involve twice as many values as the averages that created A1 . Therefore, the second averaged signal reﬂects more long term trends than those reﬂected in the ﬁrst averaged signal. We also found in Section 2.2 that this signal f = (4, 6, 10, 12, 8, 6, 5, 5) has the second ﬂuctuation d2 = (−6, 2). Consequently 1 1 −1 −1 1 1 −1 −1 2 , , , , 0, 0, 0, 0 + 2 0, 0, 0, 0, , , , D = −6 2 2 2 2 2 2 2 2 = (−3, −3, 3, 3, 1, 1, −1, −1). We found above that D1 = (−1, 1, −1, 1, 1, −1, 0, 0). Hence f = A2 + D2 + D1 = (8, 8, 8, 8, 6, 6, 6, 6) + (−3, −3, 3, 3, 1, 1, −1, −1) + (−1, 1, −1, 1, 1, −1, 0, 0). This formula further illustrates the idea of MRA. The full resolution signal f is produced from a very low resolution, averaged signal A2 consisting of repetitions of the two averaged values, 8 and 6, to which are added two detail signals. The ﬁrst addition supplements this averaged signal with enough details © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.5 Signal compression

21

to produce the next higher resolution averaged signal (5, 5, 11, 11, 7, 7, 5, 5), and the second addition then supplies enough further details to produce the full resolution signal f . In general, if the number N of signal values is divisible k times by 2, then a k-level MRA: f = Ak + Dk + · · · + D2 + D1 can be performed on the signal f . Rather than subjecting the reader to the gory details, we conclude by describing a computer example generated using FAWAV. In Figure 2.3 we show a 10-level Haar MRA of the signal f shown in Figure 2.1(a). This signal has 210 values so 10 levels of MRA are possible. On the top of Figure 2.3(a), the graph of A10 is shown; it consists of a single value repeated 210 times. This value is the average of all 210 values of the signal f . The graph directly below it is of the signal A9 which equals A10 plus the details in D10 . Each successive averaged signal is shown, from A10 through A1 . By successively adding on details, the full signal in Figure 2.1(a) is systematically constructed.

2.5

Signal compression

In Section 2.2 we saw that the Haar transform can be used to localize the energy of a signal into a shorter subsignal. In this section we show how this redistribution of energy can be used to compress audio signals. By compressing an audio signal we mean converting the signal data into a new format that requires less bits to transmit. When we use the term, audio signal, we are speaking somewhat loosely. Many of the signals we have in mind are indeed the result of taking discrete samples of a sound signal—as in the data in a computer audio ﬁle, or on a compact disc—but the techniques developed here

0

0.25

0.5

(a)

0.75

1

0

0.25

0.5

0.75

1

(b)

FIGURE 2.3 Haar MRA of the signal in Figure 2.1(a). The graphs are of the ten averaged signals A10 through A1 . Beginning with signal A10 on the top left down to A6 on the bottom left, then A5 on the top right down to A1 on the bottom right.

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

22

also apply to digital data transmissions and to other digital signals, such as digitized electrocardiograms or digitized electroencephalograms. There are two basic categories of compression techniques. The ﬁrst category is lossless compression. Lossless compression methods achieve completely error free decompression of the original signal. Typical lossless methods are Huﬀman compression, LZW compression, arithmetic compression, or run-length compression. These techniques are used in popular lossless compression programs, such as the kind that produce .zip ﬁles. Unfortunately, the compression ratios that can be obtained with lossless methods are rarely more than 2:1 for audio ﬁles consisting of music or speech. The second category is lossy compression. A lossy compression method is one which produces inaccuracies in the decompressed signal. Lossy techniques are used when these inaccuracies are so small as to be imperceptible. The advantage of lossy techniques over lossless ones is that much higher compression ratios can be attained. With wavelet compression methods, if we are willing to accept slight but imperceptible inaccuracies in the decompressed signal, then we can obtain compression ratios of 10:1, or 20:1, or as high as 50:1 or even 100:1. In order to illustrate the general principles of wavelet compression of signals, we shall examine, in a somewhat simpliﬁed way, how the Haar wavelet transform can be used to compress some test signals. For example, Signal 1 in Figure 2.4(a) can be very eﬀectively compressed using the Haar transform. Although Signal 1 is not a very representative audio signal, it is representative of a portion of a digital data transmission. This signal has 1024 values equally spaced over the time interval [0, 20). Most of these values are constant over long stretches, and that is the principal reason that Signal 1 can be compressed eﬀectively with the Haar transform. Signal 2 in Figure 2.5(a), however, will not compress nearly so well; this signal requires the more sophisticated wavelet transforms described in the next chapter. The basic steps for wavelet compression are as follows: Method of Wavelet Transform Compression Step 1. Perform a wavelet transform of the signal. Step 2. Set equal to 0 all values of the wavelet transform which are

insigniﬁcant, i.e., which lie below some threshold value. Step 3. Transmit only the signiﬁcant, non-zero values of the transform obtained from Step 2. This should be a much smaller data set than the original signal. Step 4. At the receiving end, perform the inverse wavelet transform of the data transmitted in Step 3, assigning zero values to the insigniﬁcant values which were not transmitted. This decompression step produces an approximation of the original signal. © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.5 Signal compression

0

5

10

15

4

4

2

2

0

0

−2

−2

−4 20

(a)

0

23

0

5

10

15

−4 20

(b)

5

10

(c)

15

1.01

4

1.005

2

1

0

0.995

−2

0.99 20

0

5

10

15

−4 20

(d)

FIGURE 2.4 (a) Signal 1. (b) 10-level Haar transform of Signal 1. (c) Energy map of Haar transform. (d) 20:1 compression of Signal 1, 100% of energy.

In this chapter we shall illustrate this method using the Haar wavelet transform. This initial discussion will be signiﬁcantly deepened and generalized in the next chapter when we discuss this method of compression in terms of various Daubechies wavelet transforms. Let’s now examine a Haar wavelet transform compression of Signal 1. We begin with Step 1. Since Signal 1 consists of 1024 = 210 values, we can perform 10 levels of the Haar transform. This 10-level Haar transform is shown in Figure 2.4(b). Notice how a large portion of the Haar transform’s values are 0, or very near 0, in magnitude. This fact provides the fundamental basis for performing an eﬀective compression. In order to choose a threshold value for Step 2, we proceed as follows. First, we arrange the magnitudes of the values of the Haar transform so that they are in decreasing order: L1 ≥ L2 ≥ L3 ≥ · · · ≥ LN where L1 is the largest absolute value of the Haar transform, L2 is the next © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

24

2. Haar wavelets

largest, etc. (In the event of a tie, we just leave those magnitudes in their original order.) We then compute the cumulative energy proﬁle of this new signal: 2 L1 L21 + L22 L21 + L22 + L23 , , , ...,1 . Ef Ef Ef For Signal 1, we show a graph of this energy proﬁle—which we refer to as the energy map of the Haar transform—in Figure 2.4(c). Notice that the energy map very quickly reaches its maximum value of 1. In fact, using FAWAV we ﬁnd that L21 + L22 + · · · + L251 = 0.999996. Ef Consequently, if we choose a threshold T that is less than L51 = 0.3536, then the values of the transform that survive this threshold will account for essentially 100% of the energy of Signal 1. We now turn to Step 3. In order to perform Step 3—transmitting only the signiﬁcant transform values—an additional amount of information must be sent which indicates the positions of these signiﬁcant transform values in the thresholded transform. This information is called the signiﬁcance map. The values of this signiﬁcance map are either 1 or 0: a value of 1 if the corresponding transform value survived the thresholding, a value of 0 if it did not. The signiﬁcance map is therefore a string of N bits, where N is the length of the signal. For the case of Signal 1, with a threshold of 0.35, there are only 51 non-zero bits in the signiﬁcance map out of a total of 1024 bits. Therefore, since most of this signiﬁcance map consists of long stretches of zeros, it can be very eﬀectively compressed using one of the lossless compression algorithms mentioned above. This compressed string of bits is then transmitted along with the non-zero values of the thresholded transform. Finally, we arrive at Step 4. At the receiving end, the signiﬁcance map is used to insert zeros in their proper locations in between the non-zero values in the thresholded transform, and then an inverse transform is computed to produce an approximation of the signal. For Signal 1 we show the approximation that results from using a threshold of 0.35 in Figure 2.4(d). This approximation used only 51 transform values; so it represents a compression of Signal 1 by a factor of 1024:51, i.e., a compression factor of 20:1. Since the compressed signal contains nearly 100% of the energy of the original signal, it is a very good approximation. In fact, the maximum error over all values is no more than 1.95 × 10−3 . Life would be simpler if the Haar transform could be used so eﬀectively for all signals. Unfortunately, if we try to use the Haar transform for threshold compression of Signal 2 in Figure 2.5(a), we get poor results. This signal, when played over a computer sound system, produces a sound similar to two low notes played on a clarinet. It has 4096 = 212 values; so we can perform 12 levels of the Haar transform. In Figure 2.5(b) we show a plot of the 12-level Haar transform of Signal 2. It is clear from this plot that a large fraction of © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.5 Signal compression

0

0.5

1

1.5

2

4

1

2

0

0

−1

−2

−2 2

(a)

0

25

0

0.5

1

1.5

−4 2

(b)

0.5

1

(c)

1.5

1.01

2

1.005

1

1

0

0.995

−1

0.99 2

0

0.5

1

1.5

−2 2

(d)

FIGURE 2.5 (a) Signal 2. (b) 12-level Haar transform of Signal 2. (c) Energy map of Haar transform. (d) 10:1 compression of Signal 2, 99.6% of energy.

the Haar transform values have signiﬁcant magnitude, signiﬁcant enough that they are visible in the graph. In fact, the energy map for the transform of Signal 2, shown in Figure 2.5(c), exhibits a much slower increase towards 1 in comparison with the energy map for the transform of Signal 1. Therefore, many more transform values are needed in order to capture a high percentage of the energy of Signal 2. In Figure 2.5(d), we show a 10:1 compression of Signal 2 which captures 99.6% of the energy of Signal 2. Comparing this compression with the original signal we see that it is a fairly poor approximation. Many of the signal values are clumped together in the compressed signal, producing a very ragged or jumpy approximation of the original signal. When this compressed version is played on a computer sound system, it produces a screechy “metallic” version of the two clarinet notes, which is not a very satisfying result. As a rule of thumb, we must capture at least 99.99% of the energy of the signal in order to produce an acceptable approximation, i.e., an approximation that is not perceptually diﬀerent from the original. Achieving this accurate an approximation for Signal 2 requires at least 1782 transform © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

26

values. Because Signal 2 itself has 4096 values, this is a compression ratio of only about 2.3:1, which is not very high. We shall see in the next chapter that Signal 2 can be compressed very eﬀectively, but we shall need more high powered wavelet transforms to do it.

2.5.1

A note on quantization

The most serious oversimpliﬁcation that we made in the discussion above is that we ignored the issue known as quantization. The term quantization is used whenever it is necessary to take into account the ﬁnite precision of numerical data handled by digital methods. For example, the numerical data used to generate the graphs of Signals 1 and 2 above were IEEE double precision numbers that use 8 bytes = 64 bits for each number. In order to compress this data even further, we can represent the wavelet transform coeﬃcients using less bits. We shall address this issue of quantization in the next chapter when we look again at the problem of compression.

2.6

Removing noise

In this section we shall begin our treatment of one of the most important aspects of signal processing, the removal of noise from signals. Our discussion in this section will introduce the fundamental ideas involved in the context of the Haar transform. In the next chapter we shall considerably deepen and generalize these ideas, in the context of the more powerful Daubechies wavelet transforms. When a signal is received after transmission over some distance, it is frequently contaminated by noise. The term noise refers to any undesired change that has altered the values of the original signal. The simplest model for acquisition of noise by a signal is additive noise, which has the form (contaminated signal) = (original signal) + (noise).

(2.30)

We shall represent this equation in a more compact way as f =s+n

(2.31)

where f is the contaminated signal, s is the original signal, and n is the noise. There are several kinds of noise. A few of the commonly encountered types are the following: 1. Random noise. The noise signal is highly oscillatory, its values switching rapidly between values above and below an average, or mean, value. For simplicity, we shall examine random noise with a mean value of 0. 2. Pop noise. This type of noise is heard on old analog recordings obtained from phonograph records. The noise is perceived as randomly occurring, isolated “pops.” As a model for this type of noise we add a few non-zero values to the original signal at isolated locations. © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.6 Removing noise

27

3. Localized random noise. Sometimes the noise appears as in type 1, but only over a short segment or segments of the signal. This occurs when there is a short-lived disturbance in the environment during transmission of the signal. Of course, there can also be noise signals which combine aspects of each of these types. In this section we shall examine only the ﬁrst type of noise, random noise. Other types will be considered later. Our approach will be similar to how we treated compression in the last section; we shall examine how noise removal is performed on two test signals using the Haar transform. For the ﬁrst test signal, the Haar transform is used very eﬀectively for removing the noise. For the second signal, however, the Haar transform performs poorly, and we shall need to use more sophisticated wavelet transforms to remove the noise from this signal. The essential principles, however, underlying these more sophisticated wavelet methods are the same principles we describe here for the Haar transform. Here is our basic method for removing random noise. Threshold Method of Wavelet Denoising Suppose that the contaminated signal f equals the transmitted signal s plus the noise signal n. Also suppose that these two conditions hold: 1. The energy of s is captured, to a high percentage, by transform values whose magnitudes are all greater than a threshold Ts > 0. 2. The noise signal’s transform values all have magnitudes which lie below a noise threshold Tn satisfying Tn < Ts . Then the noise in f can be removed by thresholding its transform: All values of its transform whose magnitudes lie below the noise threshold Tn are set to 0 and an inverse transform is performed, providing a good approximation of f . Let’s see how this method applies to Signal A shown in Figure 2.6(a). This signal was obtained by adding random noise, whose values oscillate between ±0.1 with a mean of zero, to Signal 1 shown in Figure 2.4(a). In this case, Signal 1 is the original signal and Signal A is the contaminated signal. As we saw in the last section, the energy of Signal 1 is captured very eﬀectively by the relatively few transform values whose magnitudes lie above a threshold of 0.35. So we set Ts equal to 0.35, and condition 1 in the Denoising Method is satisﬁed. Now as for condition 2, look at the 10-level Haar transform of Signal A shown in Figure 2.6(b). Comparing this Haar transform with the Haar transform of Signal 1 in Figure 2.4(b), it is clear that the added noise has contributed a large number of small magnitude values to the transform of Signal A, while the high-energy transform values of Signal 1 are plainly visible © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

28

(although slightly altered by the addition of noise). Therefore, we can satisfy condition 2 and eliminate the noise if we choose a noise threshold of, say, Tn = 0.2. This is indicated by the two horizontal lines shown in Figure 2.6(b); all transform values lying between ±0.2 are set to 0, producing the thresholded transform shown in Figure 2.6(c). Comparing Figure 2.6(c) with Figure 2.4(b) we see that the thresholded Haar transform of the contaminated signal is a close match to the Haar transform of the original signal. Consequently, after performing an inverse transform on this thresholded signal, we obtain a denoised signal that is a close match to the original signal. This denoised signal is shown in Figure 2.6(d), and it is clearly a good approximation to Signal 1, especially considering how much noise was originally present in Signal A.

0

5

10

15

4

4

2

2

0

0

−2

−2

−4 20

(a)

0

0

5

10

15

−4 20

(b)

5

10

(c)

15

4

4

2

2

0

0

−2

−2

−4 20

0

5

10

15

−4 20

(d)

FIGURE 2.6 (a) Signal A, 210 values. (b) 10-level Haar transform of Signal A. The two horizontal lines are at ±0.2 where 0.2 is a denoising threshold. (c) Thresholded transform. (d) Denoised signal.

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.6 Removing noise

2.6.1

29

RMS Error

The eﬀectiveness of noise removal can be quantitatively measured in the following way. The Root Mean Square Error (RMS Error) of the contaminated signal f compared with the original signal s is deﬁned to be (f1 − s1 )2 + (f2 − s2 )2 + · · · + (fN − sN )2 (RMS Error) = . (2.32) N For example, for Signal A the RMS Error between it and Signal 1 is 0.057. After denoising, the RMS Error between the denoised signal and Signal 1 is 0.011, a ﬁve-fold reduction in error. This gives quantitative evidence for the eﬀectiveness of the denoising of Signal A. Summarizing this example, we can say that the denoising was eﬀective for two reasons: (1) the transform was able to compress the energy of the original signal into a few high-energy values, and (2) the added noise was transformed into low-energy values. Consequently, the high-energy transform values from the original signal stood out clearly from the low-energy noise transform values which could then be eliminated by thresholding. Unfortunately, denoising with the Haar transform is not always so eﬀective. Consider, for example, Signal B shown in Figure 2.7(a). This signal consists of Signal 2, shown in Figure 2.5(a), with random noise added. We view Signal 2 as the original signal and Signal B as the contaminated signal. As with the ﬁrst case considered above, the random noise has values that oscillate between ±0.1 with a mean of zero. In this case, however, we saw in the last section that it takes a relatively large number of transform values to capture the energy in Signal 2. Most of these transform values are of low energy, and it takes many of them to produce a good approximation of Signal 2. When the random noise is added to Signal 2, then the Haar transform, just like in the previous case, produces many small transform values which lie below a noise threshold. This is illustrated in Figure 2.7(b) where we show the 12-level Haar transform of Signal B. As can be seen by comparing Figure 2.7(b) with Figure 2.5(b), the small transform values that come from the noise obscure most of the small magnitude transform values that result from the original signal. Consequently, when a thresholding is done to remove the noise, as indicated by the horizontal lines in Figure 2.7(b), this removes many of the transform values of the original signal which are needed for an accurate approximation. This can be veriﬁed by comparing the thresholded transform shown in Figure 2.7(c) with the original signal’s transform in Figure 2.5(b). In Figure 2.7(d) we show the denoised signal obtained by inverse transforming the thresholded signal. This denoised signal is clearly an unsatisfactory approximation of the original signal. By computing RMS Errors, we can quantify this judgment. The RMS Error between Signal B and Signal 2 is 0.057, while the RMS Error between the denoised signal and Signal 2 is 0.035. This shows that the error after denoising is almost two-thirds as great as the original error. For this second test case, we can say that the denoising was not eﬀective © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

30

0

0.5

1

1.5

2

2

1

1

0

0

−1

−1

−2 2

(a)

0

0

0.5

1

1.5

−2 2

(b)

0.5

1

1.5

(c)

2

2

1

1

0

0

−1

−1

−2 2

0

0.5

1

1.5

−2 2

(d)

FIGURE 2.7 (a) Signal B, 212 values. (b) 12-level Haar transform of Signal B. The two horizontal lines are at ±0.2 where 0.2 is the denoising threshold. (c) Thresholded transform. (d) Denoised signal.

because the transform could not compress the energy of the original signal into a few high-energy values lying above the noise threshold. We shall see in the next chapter that more sophisticated wavelet transforms can attain the desired compression and thus perform nearly as well at denoising Signal B as the Haar transform did for Signal A.

2.7

Notes and references

More material on the Haar transform and its applications can be found in [3]. Besides the lossy compression method described in this chapter, the Haar transform has also played a role in lossless image compression; see [4] or [5]. 1. R.W. Hamming. (1987). Dover, New York, NY.

Numerical Methods for Scientists and Engineers.

2. L.N. Trefethen. (1998). Maxims about numerical mathematics, computers, science, and life. SIAM News, Vol. 31, No. 1. © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.8 Examples and exercises

31

3. K.R. Rao. (1982). Fast Transforms: Algorithms, Analyses, Applications. Academic Press, New York, NY. 4. M. Rabbini and P.W. Jones. (1991). Digital Image Compression Techniques. SPIE Press, Bellingham, WA. 5. V.K. Heer and H-E Reinfelder. (1990). Comparison of reversible methods for data compression. In Medical Imaging IV, 354–365. Proceedings SPIE, No. 1233.

2.8

Examples and exercises

Note: Computer exercises, designed for FAWAV, are indicated by a superscript c. For example, problem 2.1.5 is a computer exercise. A subscript s means that a solution is provided. For instance, there is a solution provided for problem 2.1.1(a). Solutions are in Appendix B.

Section 2.1 Example 2.1.A

Find the ﬁrst level Haar transform of f = (2, 2, 2, 4, 4, 4).

Solution. The average of the ﬁrst pair of values is 2, the average of the second pair of values is 3, these √ and the average of the √ third √ pair √ of values is 4. Multiplying 1 averages by 2,√we obtain a1 = (2 2, 3 2,√4 2). To√compute √ d , we ﬁnd that 2 = −√ 2, √ and d3 = d1 =√(f1 − f2 )/ 2 = 0, and d2 = (f3 − f4 )/ 2 = −2/ √ √ (f5 − f6 )/ 2 = 0. Thus the ﬁrst level Haar transform of f is (2 2, 3 2, 4 2 | 0, − 2, 0). Example 2.1.B For the signal f = (2, 2, 2, 4, 4, 8), compute an approximate signal e f by inverse transforming the compressed Haar transform (a1 | 0, 0, 0) obtained by setting all the ﬂuctuation values equal to zero. Find the largest error between each value of f and e f. √ √ √ Solution. We ﬁnd that a1 = (2 2, 3 2, 6 2). Therefore, the inverse Haar transform produces e f = (2, 2, 3, 3, 6, 6). The largest error between each value of f and e f is 2. Example 2.1.C [Figure 2.1] To create Figure 2.1 you do the following. First, choose New 1-dim from FAWAV’s menu, and then choose Graph/Plot. Plot the formula3 20 x^2 (1-x)^4 cos(12 pi x) over the interval of type [0, L] with L = 1 using 1024 points. That produces the graph in Figure 2.1(a) (after selecting View/Display style and choosing Blank for the Grid style and Lines for the Plot style). To produce Figure 2.1(b), select Transforms/Wavelet and choose Haar as the Wavelet type with 1 entered for the Levels value. After plotting the transform, change to Lines as the plot style to get the graph shown in the ﬁgure. 2.1.1

Compute the ﬁrst trend and ﬁrst ﬂuctuation for the following signals:

(a)s

f = (2, 4, 6, 6, 4, 2)

(b)

f = (−1, 1, 2, −2, 4, −4, 2, 2)

formula is saved in the archive Primer Formulas.zip which can be downloaded from the book’s website. After extracting the formulas from this archive, you can retrieve them using the Load button under the text box in the graphing procedure. 3 This

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

32 (c)s

f = (1, 2, 3, 3, 2, 2, 1, 1)

(d)

f = (2, 2, 4, 4, 6, 6, 8, 8)

2.1.2 Given the following Haar transformed signals, ﬁnd the original signals f that correspond to them using Formula (2.6). √ √ √ √ √ √ (2 2, − 2, 3 2, − 2 | 0, 2, 0, 2) (a)s √ √ √ √ √ √ √ (b) (4 2, 3 2, − 2, 2 2 | 2, − 2, 0, 2 2) √ √ √ √ √ √ (3 2, 2 2, 2 2, 0 | 2 2, − 2, 2, 0) (c)s √ √ √ √ √ √ √ √ (d) (4 2, 5 2, 7 2, −4 2 | 2, 2 2, −2 2, 2) 2.1.3 For each of the signals f given below, compute an approximate signal e f by inverse transforming the compressed Haar transform (a1 | 0, . . . , 0) obtained by setting all the ﬂuctuation values equal to zero. In each case, ﬁnd the largest error between each value of f and e f. (a)s

f = (2, 2, 3, 3, 4, 5, 6, 6)

(b)

f = (1, 2, 3, 3, 2, 1)

(c)s

f = (2, −2, −2, −2, −2, 0, 2, 2)

(d)

f = (4, 4, 4, −1, −1, 1, 2, 2, 4, 6)

2.1.4 Consider again problem 2.1.3 above. When will there be a diﬀerence between a value of f and a value of the approximate signal e f , and when will the two signals’ values be the same? 2.1.5c Plot 1-level Haar transforms of the following functions—sampled uniformly over [0, 1] using 1024 points. (a)

f (x) = x2 (1 − x)

(b)s

f (x) = x4 (1 − x)6 cos(64πx)

(c)

(0.2 < x < 0.3) − 3(0.4 < x < 0.5) + 2(0.5 < x < 0.8)

(d)

f (x) = sgn(sin 12πx)

Section 2.2 Example 2.2.A For the signal f = (2, 2, 4, 6, 8, 10), ﬁnd the energies of its trend and ﬂuctuation subsignals and show that their sum equals the energy of f . √ √ √ √ √ Solution. The trend is a1 = (2 2, 5 2, 9 2) and ﬂuctuation is d1 = (0, − 2, − 2). The trend energy Ea1 is 8 + 50 + 162 = 220, and the ﬂuctuation energy is Ed1 = 0 + 2 + 2 = 4. Their sum is 224 and the energy of f is 4 + 4 + 16 + 36 + 64 + 100 = 224 so they are equal. Example 2.2.B Compute the percentage of compaction of the energy of f = (2, 2, 4, 6, 8, 10) by the 1-level Haar transform, in terms of Ea1 /Ef . Solution. In Example 2.2.A, we found that Ea1 = 220 and that Ef = 224. Therefore, Ea1 /Ef = 220/224 = 0.982 . . . © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.8 Examples and exercises

33

Example 2.2.C [Figure 2.2] To graph Figure 2.2(a), you plot (after selecting Edit/Points used and selecting 4096 as the number of points): 50x^2(1-x)^6cos(12pi x)(0<x

© 2008 by Taylor & Francis Group, LLC

C7451_FM.indd 1

1/3/08 11:28:56 AM

A Primer on WAVELETS and Their Scientific Applications SECOND EDITION

James S. Walker University of Wisconsin Eau Claire, Wisconsin, U.S.A.

© 2008 by Taylor & Francis Group, LLC

C7451_FM.indd 3

1/3/08 11:28:56 AM

Chapman & Hall/CRC Taylor & Francis Group 6000 Broken Sound Parkway NW, Suite 300 Boca Raton, FL 33487-2742 © 2008 by Taylor & Francis Group, LLC Chapman & Hall/CRC is an imprint of Taylor & Francis Group, an Informa business No claim to original U.S. Government works Printed in the United States of America on acid-free paper 10 9 8 7 6 5 4 3 2 1 International Standard Book Number-13: 978-1-58488-745-4 (Softcover) This book contains information obtained from authentic and highly regarded sources Reasonable efforts have been made to publish reliable data and information, but the author and publisher cannot assume responsibility for the validity of all materials or the consequences of their use. The Authors and Publishers have attempted to trace the copyright holders of all material reproduced in this publication and apologize to copyright holders if permission to publish in this form has not been obtained. If any copyright material has not been acknowledged please write and let us know so we may rectify in any future reprint Except as permitted under U.S. Copyright Law, no part of this book may be reprinted, reproduced, transmitted, or utilized in any form by any electronic, mechanical, or other means, now known or hereafter invented, including photocopying, microfilming, and recording, or in any information storage or retrieval system, without written permission from the publishers. For permission to photocopy or use material electronically from this work, please access www. copyright.com (http://www.copyright.com/) or contact the Copyright Clearance Center, Inc. (CCC) 222 Rosewood Drive, Danvers, MA 01923, 978-750-8400. CCC is a not-for-profit organization that provides licenses and registration for a variety of users. For organizations that have been granted a photocopy license by the CCC, a separate system of payment has been arranged. Trademark Notice: Product or corporate names may be trademarks or registered trademarks, and are used only for identification and explanation without intent to infringe. Library of Congress Cataloging-in-Publication Data Walker, James S. A primer on wavelets and their scientific applications / James S. Walker. -- 2nd ed. p. cm. Includes bibliographical references and index. ISBN 978-1-58488-745-4 (alk. paper) 1. Wavelets (Mathematics) I. Title. QA403.3.W33 2008 515’.2433--dc22

2007048858

Visit the Taylor & Francis Web site at http://www.taylorandfrancis.com and the CRC Press Web site at http://www.crcpress.com © 2008 by Taylor & Francis Group, LLC

C7451_FM.indd 4

1/3/08 11:28:56 AM

i

i

i

i

To Berlina and Damian

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

About the author James S. Walker received his doctorate from the University of Illinois, Chicago, in 1982. He has been a professor of Mathematics at the University of Wisconsin-Eau Claire since 1982. He has published papers on Fourier analysis, wavelet analysis, complex variables, and logic, and is the author of four books on Fourier analysis, FFTs, and wavelet analysis. He has been married to Angela Huang since 1993, and they are the proud parents of two children, Berlina and Damian, to whom this book is dedicated.

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

Preface

“Wavelet theory” is the result of a multidisciplinary eﬀort that brought together mathematicians, physicists and engineers...this connection has created a ﬂow of ideas that goes well beyond the construction of new bases or transforms. —St´ephane Mallat1

The past two decades have witnessed an explosion of activity in wavelet analysis. Thousands of research papers have been published on the theory and applications of wavelets. Wavelets provide a powerful and remarkably ﬂexible set of tools for handling fundamental problems in science and engineering. With the second edition of this primer, I aim to provide the most up-to-date introduction to this exciting ﬁeld. What is new in this second edition? I have added a lot of new material to this second edition. The most important additions are the following: • Exercises and additional examples: At the end of the chapters, I provide a total of 104 worked examples, and 222 exercises (with solutions provided for representative problems). These examples and exercises constitute a book-in-itself of review material, and should make this primer a valuable text for courses on wavelets. • Biorthogonal wavelets: I have added two sections that describe the most important examples of these wavelets, and I describe their application to image compression. 1 Mallat’s

quote is from [4] listed in the Notes and references section for Chapter 3.

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

• Image compression: I have included a mini-course on image compression, with descriptions of two state-of-the art algorithms (ASWDR and JPEG 2000). This discussion includes a tutorial on arithmetic compression, an important topic that is omitted from most (if not all) introductions to image compression. • Image denoising: I have included a lot more material on image denoising, including a description of a high-performance wavelet denoiser. In addition to the standard topic of removing Gaussian random noise, which most books limit themselves to, I also describe a technique for removing isolated, randomly positioned, clutter from images. • Gabor transforms and time-frequency analysis: I provide a concise introduction to the fundamentals of time-frequency analysis and show how it is applied to musical theory, musical synthesis, and audio denoising. An important new mathematical concept in musical theory, known as the multiresolution principle, is described here for the ﬁrst time in book form. I also provide an extensive list of references to the ﬁeld of timefrequency analysis, a topic that is given short shrift in most introductions to wavelets. • Project ideas: I provide a list, with brief descriptions, of a set of research projects on wavelets and time-frequency analysis. See Appendix A. • Book website: To keep the book as current as possible, and to provide a wealth of additional online material such as software, sound and image ﬁles, and links to other web resources on wavelets, I have created a website for the book. I describe this website in more detail later in the preface. • Enhanced list of references: I have more than doubled the list of references for the book; it now totals about 150 items. At the website for the book I provide an online version of this bibliography with a great many links for free downloading of important papers on wavelets and time-frequency analysis. What problems can be solved with wavelets? For an idea of the wide range of problems that can be solved with wavelets, here is a list of some of the problems that I shall discuss: • Audio denoising: Long distance telephone messages, for instance, often contain signiﬁcant amounts of noise. How do we remove this noise? • Signal compression: The eﬃcient transmission of large amounts of data, over the Internet for example, requires some kind of compression. How do we compress this data without losing signiﬁcant information? © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

• Object detection: What methods can we use to pick out a small image, say of an aircraft, from the midst of a larger more complicated image? • Image denoising: Images formed by microscopes and other devices are often contaminated by large amounts of unwanted clutter (referred to as noise). Can this noise be removed in order to clarify the image? • Image enhancement: When an optical microscope image is recorded, it often suﬀers from blurring. How can the appearance of the objects in these images be sharpened? • Image recognition: How do humans recognize faces? Can we teach machines to do it? • Diagnosing heart trouble: Is there a way to detect abnormal heartbeats, hidden within a complicated electrocardiogram? • Speech recognition: What factors distinguish consonants from vowels? How do humans recognize diﬀerent voices? • Musical analysis. How can mathematics help us better understand the nature of musical composition? All of these problems can be tackled using wavelets. I will show how during the course of this book. Why do we need another wavelet book? The goal of this primer is to guide the reader through the main ideas of wavelet analysis in order to facilitate a knowledgeable reading of the present research literature, especially in the applied ﬁelds of audio and image processing and biomedicine. Although there are many excellent books on the theory of wavelets, these books are focused on the construction of wavelets and their mathematical properties. Furthermore, they are all written at a graduate school level of mathematical and/or engineering expertise. There remains a real need for a basic introduction, a primer, which uses only elementary algebra and a smidgen of calculus to explain the underlying ideas behind wavelet analysis, and devotes the majority of its pages to explaining how these underlying ideas can be applied to solve signiﬁcant problems in audio and image processing and in biology and medicine. To keep the mathematics as elementary as possible, I focus on the discrete theory. It is in the continuous theory of wavelet analysis where the most diﬃcult mathematics lies; yet when this continuous theory is applied it is almost always converted into the discrete approach that I describe in this primer. Focusing on the discrete case will allow us to concentrate on the applications of wavelet analysis while at the same time keeping the mathematics under control. On the rare occasions when we need to use more advanced mathematics, I shall mark these discussions oﬀ from the main text by putting them © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

into subsections that are marked by asterisks in their titles. An eﬀort has been made to ensure that subsequent discussions do not rely on this more advanced material. Chapter summaries Chapter 1 summarizes our subject. Using examples from imaging and audio, I brieﬂy explain wavelet analysis in non-mathematical terms. In Chapter 2 I introduce the simplest wavelets, the Haar wavelets. I also introduce many of the basic concepts—wavelet transforms, energy conservation and compaction, multiresolution analysis, compression and denoising—that will be used in the remainder of the book. For this reason, I devote more pages to the theory of Haar wavelets than perhaps they deserve alone; keep in mind that this material will be ampliﬁed and generalized throughout the remainder of the book. In Chapter 3 I describe the Daubechies wavelets, which have played a key role in the explosion of activity in wavelet analysis. After a brief introduction to their mathematical properties, I describe several applications of these wavelets. First, I explain in detail how they can be used to compress audio signals—this application is vital to the ﬁelds of telephony and telecommunications. Second, I describe how the method of thresholding provides a powerful technique for removing random noise (static) from audio signals. Removing random noise is a fundamental necessity when dealing with all kinds of data in science and engineering. The threshold method, which is analogous to how our nervous system responds only to inputs above certain thresholds, provides a nearly optimal method for removing random noise. Besides random noise, Daubechies wavelets can also be used to remove isolated “pop-noise” from audio. The chapter concludes with an introduction to biorthogonal wavelets, which have played an important part in wavelet image compression algorithms such as JPEG 2000. Wavelet analysis can also be applied to images. In Chapter 4 I provide a mini-course on wavelet-based image compression. I describe two state-of-theart image compression algorithms, including the new JPEG 2000 standard. I also discuss wavelet-based image denoising, including examples motivated by optical microscopy and scanning tunnelling microscopy. Chapter 4 concludes with some examples from image processing. I examine edge detection, and the sharpening of blurred images, and an example from computer vision where wavelet methods can be used to enormously increase the speed of identiﬁcation of an image. Chapter 5 relates wavelet analysis to frequency analysis. Frequency analysis, also known as Fourier analysis, has long been one of the cornerstones of the mathematics of science and engineering. I shall brieﬂy describe how wavelets are characterized in terms of their eﬀects on the frequency content of signals. One application that I discuss is object identiﬁcation—locating a small object within a complicated scene—where wavelet analysis in concert © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

with Fourier analysis provides a powerful approach. The chapter concludes with a short introduction to the theory and application of Gabor transforms, an important relative of wavelet transforms. In recent years Gabor transforms have been increasingly used in sound processing and denoising, where they may be even more useful than wavelet transforms. In the ﬁnal chapter I deal with some extensions which reach beyond the fundamentals of wavelets. I describe a generalization of wavelet transforms known as wavelet packet transforms. I apply these wavelet packet transforms to compression of audio signals, images, and ﬁngerprints. Then I turn to the subject of continuous wavelet transforms, as they are implemented in a discrete form on a computer. Continuous wavelet transforms are widely used in seismology and have also been used very eﬀectively for analyzing speech and electrocardiograms. We shall also see how they apply to analyzing musical rhythm. Unfortunately, for reasons of space, I could not examine several important aspects of wavelet theory, such as local cosine series, best-basis algorithms, or fast matrix multiplication. At the end of the ﬁnal chapter, in its Notes and references section, I provide a guide to the literature for exploring omitted topics. Computer software Without question the best way to learn about applications of wavelets is to experiment with making such applications. This experimentation is typically done on a computer. In order to simplify this computer experimentation, I have created free software, called FAWAV. Further details about FAWAV can be found in Appendix C; suﬃce it to say for now that it is designed to allow the reader to duplicate all of the applications described in this primer and to experiment with other ideas. In this second edition I have added several new features, including image compression and denoising of both images and audio, which should add to the usefulness of FAWAV for digital signal processing. Notes and references In the Notes and references sections that conclude each chapter, I provide the reader with ample references where further information can be found. These references are numbered consecutively within each chapter. At the risk of some duplication, but with the reader’s convenience in mind, I have also compiled all references into a comprehensive Bibliography. Primer website To keep the material in the book as up-to-date as possible, and to facilitate classroom use, there is a website for this book at the following address: http://www.uwec.edu/walkerjs/Primer/ © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

This website contains a lot of supplementary material for the book, including the following: • The latest version of FAWAV and its User Manual. The User Manual is a comprehensive introduction to the use of FAWAV for digital signal processing. • Papers for downloading, including links to all papers cited in the primer. • Sound and image ﬁles for use with FAWAV, including all the sound and image ﬁles referred to in the primer. • A collection of Project ideas for carrying out your own research. This online collection will include an archive of preprints and publications that deal with these projects, and supplement the list from the primer with new topics. • Links to other websites dealing with wavelet and time-frequency analysis, including all the websites listed in the primer. Acknowledgments It is a pleasure to acknowledge several people who have helped me in writing this book. First, I must thank my wife Angela for her constant love and support. I would also like to thank my executive editor, Bob Stern, for his encouragement and help. I appreciate several suggestions from Steve Krantz that helped me put some ﬁnishing touches on the book. Thanks to Thomas Maresca for some enlightening discussion on FAWAV and its use in this primer. Thanks also to the UWEC Oﬃce of University Research and Creative Activities for a grant that gave me some extra time away from my teaching responsibilities for work on this project. Finally, my students in my Digital Signal Processing and Digital Image Processing courses deserve a warm thank you for all of their many suggestions and questions, which helped me to write a better book. I especially would like to thank Xiaowen Cheng, Jarod Hart, Nicholas Nelson, Derek Olson, and Elizabeth Otteson for their contributions. James S. Walker Eau Claire, Wisconsin

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

Contents 1 Overview 1.1 What is wavelet analysis? . . . . . . . . . . . . . . . . . . . . . 1.2 Notes and references . . . . . . . . . . . . . . . . . . . . . . . . 2 Haar wavelets 2.1 The Haar transform . . . . . . . . . . . . . . . 2.1.1 Haar transform, 1-level . . . . . . . . . 2.2 Conservation and compaction of energy . . . . 2.2.1 Conservation of energy . . . . . . . . . . 2.2.2 Haar transform, multiple levels . . . . . 2.2.3 Justiﬁcation of conservation of energy . 2.3 Haar wavelets . . . . . . . . . . . . . . . . . . . 2.4 Multiresolution analysis . . . . . . . . . . . . . 2.4.1 Multiresolution analysis, multiple levels 2.5 Signal compression . . . . . . . . . . . . . . . . 2.5.1 A note on quantization . . . . . . . . . 2.6 Removing noise . . . . . . . . . . . . . . . . . . 2.6.1 RMS Error . . . . . . . . . . . . . . . . 2.7 Notes and references . . . . . . . . . . . . . . . 2.8 Examples and exercises . . . . . . . . . . . . .

1 1 4

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

5 6 7 9 10 11 12 14 16 19 21 26 26 29 30 31

3 Daubechies wavelets 3.1 The Daub4 wavelets . . . . . . . . . . . . . . . . . . 3.1.1 Remarks on small ﬂuctuation values ∗ . . . . 3.2 Conservation and compaction of energy . . . . . . . 3.2.1 Justiﬁcation of conservation of energy ∗ . . . 3.2.2 How wavelet and scaling numbers are found ∗ 3.3 Other Daubechies wavelets . . . . . . . . . . . . . . 3.3.1 Coiﬂets . . . . . . . . . . . . . . . . . . . . . 3.4 Compression of audio signals . . . . . . . . . . . . . 3.4.1 Quantization and the signiﬁcance map . . . . 3.5 Quantization, entropy, and compression . . . . . . . 3.6 Denoising audio signals . . . . . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

41 41 49 50 50 53 54 58 61 62 65 69

. . . . . . . . . . . . . . .

. . . . . . . . . . . . . . .

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

3.6.1 Choosing a threshold value . . . . . 3.6.2 Removing pop noise and background 3.7 Biorthogonal wavelets . . . . . . . . . . . . 3.7.1 Daub 5/3 system . . . . . . . . . . . 3.7.2 Daub 5/3 inverse . . . . . . . . . . . 3.7.3 MRA for the Daub 5/3 system . . . 3.7.4 Daub 5/3 transform, multiple levels 3.7.5 Daub 5/3 integer-to-integer system . 3.8 The Daub 9/7 system . . . . . . . . . . . . 3.9 Notes and references . . . . . . . . . . . . . 3.10 Examples and exercises . . . . . . . . . . .

. . . . static . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

70 73 75 76 78 78 80 82 83 85 87

4 Two-dimensional wavelets 4.1 Two-dimensional wavelet transforms . . . . . . . . . . . . 4.1.1 Discrete images . . . . . . . . . . . . . . . . . . . . 4.1.2 2D wavelet transforms . . . . . . . . . . . . . . . . 4.1.3 2D wavelets and scaling images . . . . . . . . . . . 4.2 Compression of images—fundamentals . . . . . . . . . . . 4.2.1 Error measures in image processing . . . . . . . . . 4.2.2 Comparing JPEG with wavelet-based compressors 4.3 Fingerprint compression . . . . . . . . . . . . . . . . . . . 4.4 The WDR algorithm . . . . . . . . . . . . . . . . . . . . . 4.4.1 Bit-plane encoding . . . . . . . . . . . . . . . . . . 4.4.2 Diﬀerence reduction . . . . . . . . . . . . . . . . . 4.4.3 Arithmetic compression . . . . . . . . . . . . . . . 4.5 The ASWDR algorithm . . . . . . . . . . . . . . . . . . . 4.5.1 Arithmetic compression . . . . . . . . . . . . . . . 4.5.2 Relation to vision . . . . . . . . . . . . . . . . . . . 4.6 Important image compression features . . . . . . . . . . . 4.6.1 Progressive transmission/reconstruction . . . . . . 4.6.2 Lossless compression . . . . . . . . . . . . . . . . . 4.6.3 Region-of-interest . . . . . . . . . . . . . . . . . . . 4.7 JPEG 2000 image compression . . . . . . . . . . . . . . . 4.7.1 Compressing color images . . . . . . . . . . . . . . 4.8 Denoising images . . . . . . . . . . . . . . . . . . . . . . . 4.8.1 The TAWS algorithm . . . . . . . . . . . . . . . . 4.8.2 Comparison with Wiener denoising . . . . . . . . . 4.8.3 Estimation of noise standard deviation ∗ . . . . . . 4.8.4 Removal of clutter noise . . . . . . . . . . . . . . . 4.9 Some topics in image processing . . . . . . . . . . . . . . 4.9.1 Edge detection . . . . . . . . . . . . . . . . . . . . 4.9.2 Edge enhancement . . . . . . . . . . . . . . . . . . 4.9.3 Image recognition . . . . . . . . . . . . . . . . . . 4.10 Notes and references . . . . . . . . . . . . . . . . . . . . . 4.11 Examples and exercises . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

97 97 98 99 102 104 107 108 110 113 113 116 119 123 125 126 127 127 128 130 130 132 133 133 134 136 137 139 139 140 141 144 147

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

5 Frequency analysis 5.1 Discrete Fourier analysis . . . . . . . . . . . . . . . . . 5.1.1 Frequency content of wavelets . . . . . . . . . . 5.2 Deﬁnition of the DFT and its properties . . . . . . . . 5.2.1 Properties of the DFT . . . . . . . . . . . . . . 5.2.2 z-transforms ∗ . . . . . . . . . . . . . . . . . . . 5.3 Frequency description of wavelet analysis . . . . . . . 5.3.1 Low-pass and high-pass ﬁltering ∗ . . . . . . . . 5.4 Correlation and feature detection . . . . . . . . . . . . 5.4.1 DFT method of computing correlations . . . . 5.4.2 Proof of DFT eﬀect on correlation ∗ . . . . . . 5.4.3 Normalized correlations and feature detection ∗ 5.5 Object detection in 2D images . . . . . . . . . . . . . 5.6 Creating scaling signals and wavelets ∗ . . . . . . . . . 5.7 Gabor transforms and spectrograms . . . . . . . . . . 5.8 Musical analysis . . . . . . . . . . . . . . . . . . . . . 5.8.1 Analysis of Stravinsky’s Firebird Suite . . . . 5.8.2 Analysis of a Chinese folk song . . . . . . . . . 5.9 Inverting Gabor transforms . . . . . . . . . . . . . . . 5.10 Gabor transforms and denoising . . . . . . . . . . . . . 5.11 Notes and references . . . . . . . . . . . . . . . . . . . 5.12 Examples and exercises . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . .

167 168 169 170 171 173 174 178 180 181 183 183 185 188 192 195 197 199 201 203 206 210

6 Beyond wavelets 6.1 Wavelet packet transforms . . . . . . . . . . . . . 6.2 Wavelet packet transforms—applications . . . . . 6.2.1 Fingerprint compression . . . . . . . . . . 6.3 Continuous wavelet transforms . . . . . . . . . . 6.4 Gabor wavelets and speech analysis . . . . . . . . 6.4.1 Musical analysis: formants in song lyrics . 6.5 Percussion scalograms and musical rhythm . . . 6.5.1 Analysis of a complex percussive rhythm 6.5.2 Multiresolution Principle for rhythm . . . 6.6 Notes and references . . . . . . . . . . . . . . . . 6.6.1 Additional references . . . . . . . . . . . . 6.7 Examples and exercises . . . . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

. . . . . . . . . . . .

223 223 225 228 228 232 236 237 241 241 241 242 246

A Projects A.1 Music . . . . . . . . . . . A.2 Noise removal from audio A.3 Wavelet image processing A.4 References . . . . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

255 255 256 256 257

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

B Selected exercise B.1 Introduction . B.2 Chapter 2 . . B.3 Chapter 3 . . B.4 Chapter 4 . . B.5 Chapter 5 . . B.6 Chapter 6 . .

solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

259 259 259 262 268 273 279

C Wavelet software 283 C.1 Installing the book’s software . . . . . . . . . . . . . . . . . . . 284 C.2 Other software . . . . . . . . . . . . . . . . . . . . . . . . . . . 284 C.3 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 285 Bibliography

287

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

Chapter 1

Overview

Thought is impossible without an image. —Aristotle

We provide a brief, non-mathematical, summary of what wavelet analysis is all about. Examples are discussed from imaging and audio. In the Notes and References section we indicate where some good material on the history of wavelet analysis can be found.

1.1

What is wavelet analysis?

Wavelet analysis decomposes sounds and images into component waves of varying durations, called wavelets. These wavelets are localized vibrations of a sound signal, or localized variations of detail in an image. They can be used for a wide variety of fundamental signal processing tasks, such as compression, removing noise, or enhancing a recorded sound or image in various ways. In this section we give a couple of illustrations of what wavelet analysis involves, chosen from the ﬁelds of imaging and audio. In line with the dictum of Aristotle quoted above, we shall ﬁrst consider the ﬁeld of imaging. The heart of wavelet analysis is multiresolution analysis. Multiresolution analysis is the decomposition of a signal (such as an image) into subsignals (subimages) of diﬀerent size resolution levels. As an example, consider how the artist Ray Freer describes his process of painting: My approach to watercolor painting is very similar to that with oil or pastel. I start with a broad description of the essential masses and then © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

1. Overview

2

progress through several stages of reﬁnement until I have achieved the ﬁnal result.1

To see what Freer means, look at Figure 1.1. In Figure 1.1(a) we show an initial image, corresponding to the “broad description of the essential masses.” Figure 1.1(d) shows the details that need to be added to image (a) in order to obtain the ﬁrst reﬁnement shown in Figure 1.1(b). These details consist of regions of black where image (a) is to be darkened, and white where image (a) is to be brightened, and a dull gray background where image (a) is left unchanged. If we view these regions as a surface of water (the white parts as higher, the black parts as lower, above the undisturbed gray surface) then the white and black regions in (d) are localized waves, wavelets, that represent the changes to be made to image (a). Similarly, in image (e) we show the wavelets representing the changes to be made to image (b) to produce the ﬁnal image (c). Notice that the wavelets in (e) are ﬁner in detail, smaller scale, than the wavelets in (d). In wavelet analysis, we begin with the reverse of the process shown in Figure 1.1. We start with the ﬁnal image and determine the ﬁnest scale wavelets contained in it, and remove those to give a less detailed image. We then determine the next ﬁnest scale wavelets contained in that image and remove those, and so on, until we arrive at an image consisting of “broad masses,” a very low resolution image. This process is called the wavelet transform of the image. The wavelet transform thus provides us with the recipe for synthesizing the image as we described above. What advantages do we gain from this wavelet transform process? To see what we gain in terms of the painting example just discussed, we quote from Enrico Coen’s excellent book The Art of Genes: Why doesn’t the artist go directly to the detailed last stage, circumventing the other steps? What purpose do the previous stages serve? One problem with going too quickly into painting details is that they can be easily misplaced. If you start by painting the eyes in great detail, you might ﬁnd later on that they are not quite in the right place in relation to the other features. You would have to go back and start all over again. By ﬁrst giving the overall layout and gradually narrowing down to the details, the correct relations between features can be established gradually at each scale, from the overall location of the face, to the position of the eyes, to the details within each eye. To help them achieve this, some artists screw up their eyes when looking at the subject during the early stages of a painting, blurring their vision on purpose so as to more easily see the broad arrangements and avoid being distracted by details.2

The last sentence quoted gives us one clue as to how wavelet analysis removes random noise (random clutter) from images. If you “screw up” your eyes to 1 From 2 From

[1], page 132. [1], page 133.

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

1.1 What is wavelet analysis?

(a)

(b)

(d)

3

(c)

(e)

FIGURE 1.1 Drawing a portrait. Figure (a) is the first, rough draft. Figure (d) shows details that are added to (a) to create the second draft (b). In (d), points are blacker if the corresponding points in (a) are to be darkened. They are whiter if the corresponding points in (a) are to be brightened. And they are a dull gray if no change is needed in (a). Similarly, in (e) we show details that are needed to change (b) into the third, and final, draft (c).

blur your vision when looking at one of the noisy images in Chapter 4 [say Figure 4.17(b) on page 137, or Figure 4.18(a) on page 138], you will see that the noise is removed from the blurred image. Of course, there is much more to wavelet-based denoising than that. We will show in Chapter 4 how to add back in details to obtain sharply focused images with the noise removed. But what does this description of wavelet analysis have to do with compression of images? We can give some idea in the following way. Notice that the wavelets in Figure 1.1(d) are of a signiﬁcantly larger scale than the ones in Figure 1.1(e). This corresponds to relatively fewer brushstrokes (with a broader brush) needed to reﬁne the image in (a) to produce the image in (b). Likewise, to produce the ﬁrst draft image (a), even fewer and broader brushstrokes are needed. Because we need progressively fewer brushstrokes at these lower resolutions, the image can be compressed by specifying how to produce the lower resolution, ﬁrst draft, image and then how to successively reﬁne it (rather than trying to specify how to paint in all the ﬁne details in the original image). We will see in Chapter 4 that this, in essence, is what wavelet transform compression is doing. While 2D images are easier to see, we will begin our mathematical presentation in the next chapter with 1D signals because the mathematics is simpler © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

1. Overview

4

in 1D. The overall process, however, is similar to what we described above for images. The wavelet transform speciﬁes how to remove the ﬁnest details from a sound signal (the shortest duration wavelets) to produce a lower resolution signal. But what does “lower resolution” mean for a sound signal? It means a signal where the oscillations are of longer duration in comparison to the short duration wavelets subtracted out from the original. Steven Pinker describes this very well in terms of musical sounds: It [musical analysis] dissects the melody into essential parts and ornaments. The ornaments are stripped oﬀ and the essential parts further dissected into even more essential parts and ornaments on them. . . we sense it when we recognize variations of a piece in classical music or jazz. The skeleton of the melody is conserved while the ornaments diﬀer from variation to variation.3

The ﬁrst two sentences of this quotation are an excellent, non-mathematical, summary of what the wavelet transform is doing to 1D signals. Keep it in mind in the next two chapters when we discuss 1D wavelet transforms. We will begin with the Haar wavelet transform in Chapter 2. This transform provides the simplest way to decompose a 1D signal by stripping oﬀ details (called fluctuations) to obtain a skeleton (called a trend ). After setting the stage with Haar wavelets, we then turn to the Daubechies wavelet transforms in Chapter 3. These transforms provide a powerful approach to signal compression and signal denoising, and a deep understanding of signal synthesis.

1.2

Notes and references

For those readers interested in the history of wavelet analysis, a good place to start is the article by Burke [3], which has been expanded into a book [4]. There is also some history, of a more sophisticated mathematical nature, in the books by Meyer [5] and Daubechies [6]. 1. E. Coen. (1999). The Art of Genes. Oxford University Press, Oxford, UK. 2. S. Pinker. (1997). How the Mind Works. Norton, NY. 3. B. Burke. (1994). The mathematical microscope: waves, wavelets, and beyond. In A Positron Named Priscilla, Scientific Discovery at the Frontier, 196–235, Nat. Academy Press. Edited by M. Bartusiak. 4. B. Burke Hubbard. (1998). The World According to Wavelets, Second Edition. AK Peters, Wellesley, MA. 5. Y. Meyer. (1993). Wavelets. Algorithms and Applications. SIAM, Philadelphia, PA. 6. I. Daubechies. (1992). Ten Lectures on Wavelets. SIAM, Philadelphia, PA.

3 From

[2], page 533.

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

Chapter 2

Haar wavelets

The purpose of computing is insight, not numbers. —Richard W. Hamming

The purpose of computing is insight, not pictures. —Lloyd N. Trefethen1

A Haar wavelet is the simplest type of wavelet. In discrete form, Haar wavelets are related to a mathematical operation called the Haar transform. The Haar transform serves as a prototype for all other wavelet transforms. Studying the Haar transform in detail will provide a good foundation for understanding the more sophisticated wavelet transforms which we shall describe in the next chapter. In this chapter we shall describe how the Haar transform can be used for compressing audio signals and for removing noise. Our discussion of these applications will set the stage for the more powerful wavelet transforms to come and their applications to these same problems. One distinctive feature that the Haar transform enjoys is that it lends itself easily to simple hand calculations. We shall illustrate many concepts by both simple hand calculations and more involved computer computations.

1 Hamming’s

quote is from [1]. Trefethen’s quote is from [2].

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

6

2.1

The Haar transform

In this section we shall introduce the basic notions connected with the Haar transform, which we shall examine in more detail in later sections. First, we need to deﬁne the type of signals that we shall be analyzing. We shall be working extensively with discrete signals. A discrete signal is a function of time with values occurring at discrete instants. Generally we shall express a discrete signal in the form f = (f1 , f2 , . . . , fN ), where N is a positive even integer2 which we shall refer to as the length of f . The values of f are the N real numbers f1 , f2 , . . . , fN . These values are typically measured values of an analog signal g, measured at the time values t = t1 , t2 , . . . , tN . That is, the values of f are f1 = g(t1 ), f2 = g(t2 ), . . . , fN = g (tN ) .

(2.1)

For simplicity, we shall assume that the increment of time that separates each pair of successive time values is always the same. We shall use the phrase equally spaced sample values, or just sample values, when the discrete signal has its values deﬁned in this way. An important example of sample values is the set of data values stored in a computer audio ﬁle, such as a .wav ﬁle. Another example is the sound intensity values recorded on a compact disc. A non-audio example is a digitized electrocardiogram. Like all wavelet transforms, the Haar transform decomposes a discrete signal into two subsignals of half its length. One subsignal is a running average or trend; the other subsignal is a running diﬀerence or ﬂuctuation. Let’s begin with the trend. The ﬁrst trend, a1 = (a1 , a2 , . . . , aN/2 ), for the signal f is computed by taking a running average in the following way. Its ﬁrst average of the ﬁrst pair of value, a1 , is computed by taking the √ √ values of f : (f1 +f2 )/2; and then multiplying it by 2. Thus, a1 = (f1 +f2 )/ 2. Similarly, of the next pair of values its next value a2 is computed by taking the average √ √ of f : (f3 + f4 )/2; and then multiplying it by 2. Hence, a2 = (f3 + f4 )/ 2. Continuing in this way, all of the values of a1 are produced by taking averages √ of successive pairs of values of f , and then multiplying these averages by 2. A formula for the values of a1 is am =

f2m−1 + f2m √ , 2

(2.2)

for m = 1, 2, 3, . . . , N/2. For example, suppose f is deﬁned by eight values, say f = (4, 6, 10, 12, 8, 6, 5, 5). 2 If

N is odd, then a zero can be appended to the signal to make it even.

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.1 The Haar transform

7

The average of its ﬁrst two values is 5, the average of the next two values is 11, the average of the next two values is√7, and the average of the last two values is 5. Multiplying these averages by 2 , we obtain the ﬁrst trend subsignal √ √ √ √ a1 = (5 2, 11 2, 7 2, 5 2). √ You might ask: Why perform the extra step of multiplying by 2 ? Why not just take averages? These questions √ will be answered in the next section, when we show that multiplication by 2 is needed in order to ensure that the Haar transform preserves the energy of a signal. The other subsignal is called the ﬁrst ﬂuctuation. The ﬁrst ﬂuctuation of the signal f , which is denoted by d1 = (d1 , d2 , . . . , dN/2 ), is computed by taking a running diﬀerence in the following way. Its ﬁrst value, d1 , is found by taking half the√diﬀerence of the ﬁrst pair√of values of f : (f1 − f2 )/2; and multiplying it by 2. That is, d1 = (f1 −f2 )/ 2. Likewise, its next value d2 is found by taking half the √ f : (f3 − f4 )/2; √diﬀerence of the next pair of values of and multiplying it by 2. In other words, d2 = (f3 − f4 )/ 2. Continuing in this way, all of the values of d1 are produced according to the following formula: f2m−1 − f2m √ , (2.3) dm = 2 for m = 1, 2, 3, . . . , N/2. For example, for the signal considered above, f = (4, 6, 10, 12, 8, 6, 5, 5), its ﬁrst ﬂuctuation is

2.1.1

√ √ √ d1 = (− 2, − 2, 2, 0).

Haar transform, 1-level

The Haar transform is performed in several stages, or levels. The ﬁrst level is the mapping H1 deﬁned by 1 (a1 | d1 ) f −→

H

(2.4)

from a discrete signal f to its ﬁrst trend a1 and ﬁrst ﬂuctuation d1 . For example, we showed above that √ √ √ √ √ √ √ H1 (4, 6, 10, 12, 8, 6, 5, 5) −→ (5 2, 11 2, 7 2, 5 2 | − 2, − 2, 2, 0).

(2.5)

The mapping H1 in (2.4) has an inverse. Its inverse maps the transform signal (a1 | d1 ) back to the signal f , via the following formula: aN/2 + dN/2 aN/2 − dN/2 a1 + d1 a1 − d1 √ √ √ , √ ,..., , . (2.6) f= 2 2 2 2 © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

8

√ √ √ In other words,√f1 = (a1 + d1 )/ 2, f2 = (a1 − d1 )/ 2, f3 = (a2 + d2 )/ 2, f4 = (a2 − d2 )/ 2, and so on. The reader may wish to check that the formulas in (2.6) when applied to √ √ √ √ √ √ √ (5 2, 11 2, 7 2, 5 2 | − 2, − 2, 2, 0) produce (4, 6, 10, 12, 8, 6, 5, 5) thereby inverting the transform in (2.5). Let’s now consider what advantages accrue from performing the Haar transformation. These advantages will be described in more detail later in this chapter, but some basic notions can be introduced now. All of these advantages stem from the following cardinal feature of the Haar transform (a feature that will be even more prominent for the Daubechies transforms described in the next chapter): Small Fluctuations Feature. The magnitudes of the values of the ﬂuctuation subsignal are often signiﬁcantly smaller than the magnitudes of the values of the original signal. For instance, for the signal f = (4, 6, 10, 12, 8, 6, 5, 5) considered above, its eight values have an √ average of 7. On the other hand, for its ﬁrst √ magnitude √ 1 ﬂuctuation d = (− 2, − 2, 2, 0), the average of its four magnitudes is √ 0.75 2. In this case, the magnitudes of the ﬂuctuation’s values are an average of 6.6 times smaller than the magnitudes of the original signal’s values. For a second example, consider the signal shown in Figure 2.1(a). This signal was generated from 1024 sample values of the function g(x) = 20x2 (1 − x)4 cos 12πx over the interval [0, 1). In Figure 2.1(b) we show a graph of the 1-level Haar transform of this signal. The trend subsignal is graphed on the left half, over the interval [0, 0.5), and the ﬂuctuation subsignal is graphed on the right half, over the interval [0.5, 1). It is clear that a large percentage of the ﬂuctuation’s values are close to 0 in magnitude, another instance of the Small Fluctuations Feature. Notice also that the trend subsignal looks like the original signal, √ although shrunk by half in length and expanded by a factor of 2 vertically. The reason that the Small Fluctuations Feature is generally true is that typically we are dealing with signals sampled from a continuous analog signal g with a very short time increment between the samples. In other words, the equations in (2.1) hold with a small value of the time increment h = tk+1 − tk for each k = 1, 2, . . . , N − 1. If the time increment is small enough, then successive values f2m−1 = g (t2m−1 ) and f2m = g (t2m ) of the signal f will be close to each other due to the continuity of g. Consequently, the ﬂuctuation values for the Haar transform satisfy dm =

g (t2m−1 ) − g (t2m ) √ ≈ 0. 2

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.2 Conservation and compaction of energy

0

0.25

0.5

0.75

9

1

1

0.5

0.5

0

0

−0.5

−0.5

−1 1

(a)

0

0.25

0.5

0.75

−1 1

(b)

FIGURE 2.1 (a) Signal, (b) Haar transform, 1-level.

This explains why the Small Fluctuations Feature is generally true for the Haar transform. A similar analysis shows why the trend subsignal has a graph that is similar in appearance to the ﬁrst trend. If g is continuous and the time increment is very small, then g(t2m−1 ) and g(t2m ) will be close to each other. Expressing this fact as an approximation, g(t2m−1 ) ≈ g(t2m ), we obtain the following approximation for each value am of the trend subsignal √ am ≈ 2 g(t2m ). This equation shows that a1 is approximately the same as sample values of √ 2 g(x) for x = t2 , t4 , . . . , tN . In other words, it shows that the graph of the ﬁrst trend subsignal is similar in appearance to the graph of g, as we pointed out above in regard to the signal in Figure 2.1(a). We shall examine these points in more detail in the next chapter. One of the reasons that the Small Fluctuations Feature is important is that it has applications to signal compression. By compressing a signal we mean transmitting its values, or approximations of its values, by using a smaller number of bits. For example, if we were only to transmit the trend subsignal for the signal shown in Figure 2.1(a) and then perform Haar transform inversion (treating the ﬂuctuation’s values as all zeros), then we would obtain an approximation of the original signal. Since the length of the trend subsignal is half the length of the original signal, this would achieve 50% compression. We shall discuss compression in more detail in Section 2.5. Once we have performed a 1-level Haar transform, then it is easy to perform multiple-level Haar transforms. We shall discuss this in the next section.

2.2

Conservation and compaction of energy

In the previous section we deﬁned the 1-level Haar transform. In this section we shall discuss its two most important properties: (1) It conserves the ener© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

10

gies of signals; (2) It performs a compaction of the energy of signals. We shall also complete our deﬁnition of the Haar transform by showing how to extend its deﬁnition to multiple levels.

2.2.1

Conservation of energy

An important property of the Haar transform is that it conserves the energies of signals. By the energy of a signal f we mean the sum of the squares of its values. That is, the energy Ef of a signal f is deﬁned by 2 Ef = f12 + f22 + · · · + fN .

(2.7)

We shall provide some explanation for why we use the term energy for this quantity Ef in a moment. First, however, let’s look at an example of calculating energy. Suppose f = (4, 6, 10, 12, 8, 6, 5, 5) is the signal considered in Section 2.1. Then Ef is calculated as follows: Ef = 42 + 62 + · · · + 52 = 446. Furthermore, using the values for its 1-level Haar transform √ √ √ √ √ √ √ (a1 | d1 ) = (5 2, 11 2, 7 2, 5 2 | − 2, − 2, 2, 0), we ﬁnd that E(a1 | d1 ) = 25 · 2 + 121 · 2 + · · · + 2 + 0 = 446 as well. Thus the 1-level Haar transform has kept the energy constant. In fact, this is true in general: Conservation of energy. The 1-level Haar transform conserves energy. In other words, E(a1 | d1 ) = Ef for every signal f . We will explain why this conservation of energy property is true for all signals at the end of this section. Before we go any further, we should say something about why we use the term energy for the quantity Ef . The reason is that sums of squares frequently appear in physics when various types of energy are calculated. For instance, if a particle of mass m has a velocity of v = (v1 , v2 , v3 ), then its kinetic energy is (m/2)(v12 + v22 + v32 ). Hence its kinetic energy is proportional to v12 + v22 + v32 = Ev . Ignoring the constant of proportionality, m/2, we obtain the quantity Ev which we call the energy of v. While conservation of energy is certainly an important property, it is even more important to consider how the Haar transform redistributes the energy in a signal by compressing most of the energy into the trend subsignal. For example, for the signal √ f = (4, √ 6, 10, √ 12,√8, 6, 5, 5) we found in Section 2.1 that its trend a1 equals (5 2, 11 2, 7 2, 5 2). Therefore, the energy of a1 is Ea1 = 25 · 2 + 121 · 2 + 49 · 2 + 25 · 2 = 440. © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.2 Conservation and compaction of energy

11

√ √ √ On the other hand, the ﬂuctuation d1 is (− 2, − 2, 2, 0), which has energy Ed1 = 2 + 2 + 2 + 0 = 6. Thus the energy of the trend a1 accounts for 440/446 = 98.7% of the total energy of the signal. In other words, the 1-level Haar transform has redistributed the energy of f so that over 98% is concentrated into the subsignal a1 which is half the length of f . For obvious reasons, this is called compaction of energy. As another example, consider the signal f graphed in Figure 2.1(a) and its 1-level Haar transform shown in Figure 2.1(b). In this case, we ﬁnd that the energy of the signal f is 127.308 while the energy of its ﬁrst trend a1 is 127.305. Thus 99.998% of the total energy is compacted into the half-length subsignal a1 . By examining the graph in Figure 2.1(b) it is easy to see why such a phenomenal energy compaction has occurred; the values of the ﬂuctuation d1 are so small, relative to the much larger values of the trend a1 , that its energy Ed1 contributes only a small fraction of the total energy Ea1 + Ed1 . These two examples illustrate the following general principle: Compaction of energy. The energy of the trend subsignal a1 accounts for a large percentage of the energy of the transformed signal (a1 | d1 ). Compaction of energy will occur whenever the magnitudes of the ﬂuctuation’s values are signiﬁcantly smaller than the trend’s values (recall the Small Fluctuations Feature from the last section). In Section 2.5, we shall describe how compaction of energy provides a framework for applying the Haar transform to compress signals. We now turn to a discussion of how the Haar transform can be extended to multiple levels, thereby increasing the energy compaction of signals.

2.2.2

Haar transform, multiple levels

Once we have performed a 1-level Haar transform, then it is easy to repeat the process and perform multiple level Haar transforms. After performing a 1-level Haar transform of a signal f we obtain a ﬁrst trend a1 and a ﬁrst ﬂuctuation d1 . The second level of a Haar transform is then performed by computing a second trend a2 and a second ﬂuctuation d2 for the ﬁrst trend a1 only. For example, if f = (4, 6, 10, 12, 8, 6, 5, 5)√is the above, √ signal √ considered √ then we found that its ﬁrst trend is a1 = (5 2, 11 2, 7 2, 5 2). To get the 1 second trend a2 we apply Formula (2.2)√to the√values √ of√ a . That is, we 1 add successive pairs of values of a = (5 2, 11 2, 7 2, 5 2) and divide by √ 2 obtaining a2 = (16, 12). To get √ the second ﬂuctuation d2 , we subtract √ √ √ √ 1 successive pairs of values of a = (5 2, 11 2, 7 2, 5 2) and divide by 2 obtaining a2 = (−6, 2). Thus the 2-level Haar transform of f is the signal √ √ √ (a2 | d2 | d1 ) = (16, 12 | −6, 2 | − 2, − 2, 2, 0). © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

12

For this signal f , a 3-level Haar transform can also be done, and the result is √ √ √ √ √ (a3 | d3 | d2 | d1 ) = (14 2 | 2 2 | −6, 2 | − 2, − 2, 2, 0). It is interesting to calculate the energy compaction that has occurred with the 2-level and 3-level Haar transforms that we just computed. First, we know that E(a2 | d2 | d1 ) = 446 because of conservation of energy. Second, we compute that Ea2 = 400. Thus the 2-level Haar transformed signal (a2 | d2 | d1 ) has almost 90% of the total energy of f contained in the second trend a2 which is 1/4 of the length of f . This is a further compaction, or localization, of the energy of f . Furthermore, Ea3 = 392; thus a3 contains 87.89% of the total energy of f . This is even further compaction; the 3-level Haar transform (a3 | d3 | d2 | d1 ) has almost 88% of the total energy of f contained in the third trend a3 which is 1/8 the length of f . For those readers who are familiar with Quantum Theory, there is an interesting phenomenon here that is worth noting. By Heisenberg’s Uncertainty Principle, it is impossible to localize a ﬁxed amount of energy into an arbitrarily small time interval. This provides an explanation for why the energy percentage dropped from 98% to 90% when the second-level Haar transform was computed, and from 90% to 88% when the third-level Haar transform was computed. When we attempt to squeeze the energy into ever smaller time intervals, it is inevitable that some energy leaks out. As another example of how the Haar transform redistributes and localizes the energy in a signal, consider the graphs shown in Figure 2.2. In Figure 2.2(a) we show a signal, and in Figure 2.2(b) we show the 2-level Haar transform of this signal. In Figures 2.2(c) and (d) we show the respective cumulative energy proﬁles of these two signals. By the cumulative energy proﬁle of a signal f we mean the signal deﬁned by

f12 f12 + f22 f12 + f22 + f32 , , , ...,1 . Ef Ef Ef

The cumulative energy proﬁle of f thus provides a summary of the accumulation of energy in the signal as time proceeds. As can be seen from comparing the two proﬁles in Figure 2.2, the 2-level Haar transform has redistributed and localized the energy of the original signal.

2.2.3

Justification of conservation of energy

We close this section with a brief justiﬁcation of the conservation of energy property of the Haar transform. First, we observe that the terms a21 and d21 © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.2 Conservation and compaction of energy

0

0.5

1

1.5

1.5

1.5

0.75

0.75

0

0

−0.75

−0.75

−1.5 2

(a)

0

13

0

0.5

1

1.5

−1.5 2

(b)

0.5

1

1.5

(c)

1.5

1.5

1

1

0.5

0.5

0

0

−0.5 2

0

0.5

1

1.5

−0.5 2

(d)

FIGURE 2.2 (a) Signal. (b) 2-level Haar transform of signal. (c) Cumulative energy profile of signal. (d) Cumulative energy profile of 2-level Haar transform.

in the formula E(a1 | d1 ) = a21 + · · · + a2N/2 + d21 + · · · + d2N/2 add up as follows: 2 2 f1 + f2 f1 − f2 √ √ + 2 2 2 2 f 2 − 2f1 f2 + f22 f + 2f1 f2 + f2 + 1 = 1 2 2 = f12 + f22 .

a21 + d21 =

2 2 + f2m for all other values of m. Therefore, by Similarly, a2m + d2m = f2m−1 adding a2m and d2m successively for each m, we ﬁnd that 2 . a21 + · · · + a2N/2 + d21 + · · · + d2N/2 = f12 + · · · + fN

In other words, E(a1 | d1 ) = Ef , which justiﬁes the conservation of energy property. © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

14

2.3

Haar wavelets

In this section we discuss the simplest wavelets, the Haar wavelets. This material will set the stage for the more sophisticated Daubechies wavelets described in the next chapter. We begin by discussing the 1-level Haar wavelets. These wavelets are deﬁned as 1 −1 1 W1 = √ , √ , 0, 0, . . . , 0 2 2 1 −1 1 W2 = 0, 0, √ , √ , 0, 0, . . . , 0 2 2 −1 1 1 √ √ W3 = 0, 0, 0, 0 , , 0, 0, . . . , 0 2 2 .. . 1 −1 1 WN/2 = 0, 0, . . . , 0, √ , √ . (2.8) 2 2 These 1-level Haar wavelets have a number of interesting properties. First, they each have an energy of 1. Second, √ they each consist of a rapid ﬂuctuation between just two non-zero values, ±1/ 2, with an average value of 0. Hence the name wavelet. Finally, they all are very similar to each other in that they are each a translation in time by an even number of time-units of the ﬁrst Haar wavelet W11 . The second Haar wavelet W21 is a translation forward in time by two units of W11 , and W31 is a translation forward in time by four units of W11 , and so on. The reason for introducing the 1-level Haar wavelets is that we can express the 1-level ﬂuctuation subsignal in a simpler form by making use of scalar products with these wavelets. The scalar product is a fundamental operation on two signals, and is deﬁned as follows. Scalar product: The scalar product f · g of the signals f = (f1 , f2 , . . . , fN ) and g = (g1 , g2 , . . . , gN ) is deﬁned by f · g = f1 g1 + f2 g2 + · · · + fN gN .

(2.9)

Using the 1-level Haar wavelets, we can express the values for the ﬁrst ﬂuctuation subsignal d1 as scalar products. For example, d1 =

f1 − f2 √ = f · W11 . 2

Similarly, d2 = f · W21 , and so on. We can summarize Formula (2.3) in terms of scalar products with the 1-level Haar wavelets: 1 dm = f · Wm

(2.10)

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.3 Haar wavelets

15

for m = 1, 2, . . . , N/2. We can also use the idea of scalar products to restate the Small Fluctuations Feature from Section 2.1 in a more precise form. If we say that the support of each Haar wavelet is the set of two time-indices where the wavelet is non-zero, then we have the following more precise version of the Small Fluctuations Feature: Property 1. If a signal f is (approximately) constant over the support of a 1-level Haar wavelet Wk1 , then the ﬂuctuation value dk = f · Wk1 is (approximately) zero. This property will be considerably strengthened in the next chapter. Note: From now on, we shall refer to the set of time-indices m where fm = 0 as the support of a signal f . For example, the support of the signal f = (0, 0, 5, 7, −2, 0, 2, 0) consists of the indices 3, 4, 5, and 7; while the support of the Haar wavelet W21 consists of the indices 3, 4. We can also express the 1-level trend values as scalar products with certain elementary signals. These elementary signals are called 1-level Haar scaling signals, and they are deﬁned as 1 1 = √ , √ , 0, 0, . . . , 0 2 2 1 1 1 V2 = 0, 0, √ , √ , 0, 0, . . . , 0 2 2 .. . 1 1 1 = 0, 0, . . . , 0, √ , √ . VN/2 2 2 V11

(2.11)

Using these Haar scaling signals, the values a1 , . . . , aN/2 for the ﬁrst trend are expressed as scalar products: 1 am = f · Vm

(2.12)

for m = 1, 2, . . . , N/2. The Haar scaling signals are quite similar to the Haar wavelets. They all have energy 1 and have a support consisting of just two consecutive timeindices. In fact, they are all translates by an even multiple of time-units of the ﬁrst scaling signal V11 . Unlike the Haar wavelets, however, the average values of the Haar √scaling signals are not zero. In fact, they each have an average value of 1/ 2. The ideas discussed above extend to every level. For simplicity, we restrict our discussion to the second level. The 2-level Haar scaling signals are deﬁned © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

16 by

1 1 1 1 , , , , 0, 0 . . . , 0 2 2 2 2 1 1 1 1 V22 = 0, 0, 0, 0, , , , , 0, 0, . . . , 0 2 2 2 2 .. . 1 1 1 1 2 . VN/4 = 0, 0, . . . , 0, , , , 2 2 2 2 V12 =

(2.13)

These scaling signals are all translations by multiples of four time-units of the ﬁrst scaling signal V12 , and they all have energy 1 and average value 1/2. Furthermore, the values of the 2-level trend a2 are scalar products of these scaling signals with the signal f . That is, a2 satisﬁes 2 . (2.14) a2 = f · V12 , f · V22 , . . . , f · VN/4 Likewise, the 2-level Haar wavelets are deﬁned by 1 1 −1 −1 W12 = , , , , 0, 0 . . . , 0 2 2 2 2 1 1 −1 −1 W22 = 0, 0, 0, 0, , , , , 0, 0, . . . , 0 2 2 2 2 .. . 1 1 −1 −1 2 , . = 0, 0, . . . , 0, , , WN/4 2 2 2 2

(2.15)

These wavelets all have supports of length 4, since they are all translations by multiples of four time-units of the ﬁrst wavelet W12 . They also all have energy 1 and average value 0. Using scalar products, the 2-level ﬂuctuation d2 satisﬁes 2 . (2.16) d2 = f · W12 , f · W22 , . . . , f · WN/4

2.4

Multiresolution analysis

In the previous section we discussed how the Haar transform can be described using scalar products with scaling signals and wavelets. In this section we discuss how the inverse Haar transform can also be described in terms of these same elementary signals. This discussion will show how discrete signals are synthesized by beginning with a very low resolution signal and successively adding on details to create higher resolution versions, ending with a complete synthesis of the signal at the ﬁnest resolution. This is known as multiresolution analysis (MRA). MRA is the heart of wavelet analysis. © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.4 Multiresolution analysis

17

In order to make these ideas precise, we must ﬁrst discuss some elementary operations that can be performed on signals. Given two signals f = (f1 , f2 , . . . , fN ) and g = (g1 , g2 , . . . , gN ), we can perform the following algebraic operations: Addition and Subtraction: The sum f + g of the signals f and g is deﬁned by adding their values: f + g = (f1 + g1 , f2 + g2 , . . . , fN + gN ).

(2.17)

Their diﬀerence f − g is deﬁned by subtracting their values: f − g = (f1 − g1 , f2 − g2 , . . . , fN − gN ).

(2.18)

Constant multiple: A signal f is multiplied by a constant c by multiplying each of its values by c. That is, c f = (cf1 , cf2 , . . . , cfN ).

(2.19)

For example, by repeatedly applying the addition operation, we can express a signal f = (f1 , f2 , . . . , fN ) as follows: f = (f1 , 0, 0, . . . , 0) + (0, f2 , 0, 0, . . . , 0) + · · · + (0, 0, . . . , 0, fN ). Then, by applying the constant multiple operation to each of the signals on the right side of this last equation, we obtain f = f1 (1, 0, 0, . . . , 0) + f2 (0, 1, 0, 0, . . . , 0) + · · · + fN (0, 0, . . . , 0, 1). This formula is a very natural one; it amounts to expressing f as a sum of its individual values at each discrete instant of time. 0 as If we deﬁne the elementary signals V10 , V20 , . . . , VN V10 = (1, 0, 0, . . . , 0) V20 = (0, 1, 0, 0, . . . , 0) .. . 0 VN = (0, 0, . . . , 0, 1)

(2.20)

then the last formula for f can be rewritten as 0 . f = f1 V10 + f2 V20 + · · · + fN VN

(2.21)

Formula (2.21) is called the natural expansion of a signal f in terms of the 0 . We shall now show that the Haar natural basis of signals V10 , V20 , . . . , VN MRA involves expressing f as a sum of constant multiples of a diﬀerent basis set of elementary signals: the Haar wavelets and scaling signals deﬁned in the previous section. © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

18

In the previous section, we showed how to express the 1-level Haar transform in terms of wavelets and scaling signals. It is also possible to express the inverse of the 1-level Haar transform in terms of these same elementary signals. This leads to the ﬁrst level of the Haar MRA. To deﬁne this ﬁrst level Haar MRA we make use of (2.6) to express a signal f as aN/2 aN/2 a1 a1 a2 a2 f = √ , √ , √ , √ ,..., √ , √ 2 2 2 2 2 2 dN/2 −dN/2 d1 −d1 d2 −d2 . + √ , √ , √ , √ ,..., √ , √ 2 2 2 2 2 2 This formula shows that the signal f can be expressed as the sum of two signals that we shall call the ﬁrst averaged signal and the ﬁrst detail signal. That is, we have f = A1 + D1 (2.22) where the signal A1 is called the ﬁrst averaged signal and is deﬁned by aN/2 aN/2 a1 a1 a2 a2 1 A = √ , √ , √ , √ ,..., √ , √ (2.23) 2 2 2 2 2 2 and the signal D1 is called the ﬁrst detail signal and is deﬁned by dN/2 −dN/2 d1 −d1 d2 −d2 1 D = √ , √ , √ , √ ,..., √ , √ . 2 2 2 2 2 2

(2.24)

Using Haar scaling signals and wavelets, and using the basic algebraic operations with signals, the averaged and detail signals are expressible as 1 , A1 = a1 V11 + a2 V21 + · · · + aN/2 VN/2

(2.25a)

1 D1 = d1 W11 + d2 W21 + · · · + dN/2 WN/2 .

(2.25b)

Applying the scalar product formulas for the coeﬃcients in Equations (2.10) and (2.12), we can rewrite these last two formulas as follows 1 1 )VN/2 , A1 = (f · V11 )V11 + (f · V21 )V21 + · · · + (f · VN/2 1 1 )WN/2 . D1 = (f · W11 )W11 + (f · W21 )W21 + · · · + (f · WN/2

These formulas show that the averaged signal is a combination of Haar scaling signals, with the values of the ﬁrst trend subsignal as coeﬃcients; and that the detail signal is a combination of Haar wavelets, with the values of the ﬁrst ﬂuctuation subsignal as coeﬃcients. As an example of these ideas, consider the signal f = (4, 6, 10, 12, 8, 6, 5, 5). © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.4 Multiresolution analysis

19

In Section 2.1 we found that its ﬁrst trend subsignal was √ √ √ √ a1 = (5 2, 11 2, 7 2, 5 2). Applying Formula (2.23), the averaged signal is A1 = (5, 5, 11, 11, 7, 7, 5, 5).

(2.26)

Notice how the ﬁrst averaged signal consists of the repeated average values 5, 5, and 11, 11, and 7, 7, and 5, 5 about which the values of f ﬂuctuate. Using Formula (2.25a), the ﬁrst averaged signal can also be expressed in terms of Haar scaling signals as √ √ √ √ A1 = 5 2 V11 + 11 2 V21 + 7 2 V31 + 5 2 V41 . Comparing these last two equations we can see that the positions of the repeated averages correspond precisely with the supports of the scaling signals. We also √ found √ √in Section 2.1 that the ﬁrst ﬂuctuation signal for f was d1 = (− 2, − 2, 2, 0). Formula (2.24) then yields D1 = (−1, 1, −1, 1, 1, −1, 0, 0). Thus, using the result for A1 computed above, we have f = (5, 5, 11, 11, 7, 7, 5, 5) + (−1, 1, −1, 1, 1, −1, 0, 0). This equation illustrates the basic idea of MRA. The signal f is expressed as a sum of a lower resolution, or averaged, signal (5, 5, 11, 11, 7, 7, 5, 5) added with a signal (−1, 1, −1, 1, 1, −1, 0, 0) made up of ﬂuctuations or details. These ﬂuctuations provide the added details necessary to produce the full resolution signal f . For this example, using Formula (2.25b), the ﬁrst detail signal can also be expressed in terms of Haar wavelets as √ √ √ D1 = − 2 W11 − 2 W21 + 2 W31 + 0 W41 . This formula shows that the values of D1 occur in successive pairs of rapidly ﬂuctuating values positioned at the supports of the Haar wavelets.

2.4.1

Multiresolution analysis, multiple levels

In the discussion above, we described the ﬁrst level of the Haar MRA of a signal. This idea can be extended to further levels, as many levels as the number of times that the signal length can be divided by 2. The second level of a MRA of a signal f involves expressing f as f = A 2 + D2 + D1 .

(2.27)

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

20

Here A2 is the second averaged signal and D2 is the second detail signal. Comparing Formulas (2.22) and (2.27) we see that A1 = A2 + D2 .

(2.28)

This formula expresses the fact that computing the second averaged signal A2 and second detail signal D2 simply consists of performing a ﬁrst level MRA of the signal A1 . Because of this, it follows that the second level averaged signal A2 satisﬁes 2 2 )VN/4 A2 = (f · V12 )V12 + (f · V22 )V22 + · · · + (f · VN/4

and the second level detail signal D2 satisﬁes 2 2 )WN/4 . D2 = (f · W12 )W12 + (f · W22 )W22 + · · · + (f · WN/4

For example, if f = (4, 6, 10, 12, 8, 6, 5, 5), then we found in Section 2.2 that a2 = (16, 12). Therefore 1 1 1 1 1 1 1 1 , , , , 0, 0, 0, 0 + 12 0, 0, 0, 0, , , , A2 = 16 2 2 2 2 2 2 2 2 = (8, 8, 8, 8, 6, 6, 6, 6).

(2.29)

It is interesting to compare the equations in (2.26) and (2.29). The second averaged signal A2 has values created from averages that involve twice as many values as the averages that created A1 . Therefore, the second averaged signal reﬂects more long term trends than those reﬂected in the ﬁrst averaged signal. We also found in Section 2.2 that this signal f = (4, 6, 10, 12, 8, 6, 5, 5) has the second ﬂuctuation d2 = (−6, 2). Consequently 1 1 −1 −1 1 1 −1 −1 2 , , , , 0, 0, 0, 0 + 2 0, 0, 0, 0, , , , D = −6 2 2 2 2 2 2 2 2 = (−3, −3, 3, 3, 1, 1, −1, −1). We found above that D1 = (−1, 1, −1, 1, 1, −1, 0, 0). Hence f = A2 + D2 + D1 = (8, 8, 8, 8, 6, 6, 6, 6) + (−3, −3, 3, 3, 1, 1, −1, −1) + (−1, 1, −1, 1, 1, −1, 0, 0). This formula further illustrates the idea of MRA. The full resolution signal f is produced from a very low resolution, averaged signal A2 consisting of repetitions of the two averaged values, 8 and 6, to which are added two detail signals. The ﬁrst addition supplements this averaged signal with enough details © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.5 Signal compression

21

to produce the next higher resolution averaged signal (5, 5, 11, 11, 7, 7, 5, 5), and the second addition then supplies enough further details to produce the full resolution signal f . In general, if the number N of signal values is divisible k times by 2, then a k-level MRA: f = Ak + Dk + · · · + D2 + D1 can be performed on the signal f . Rather than subjecting the reader to the gory details, we conclude by describing a computer example generated using FAWAV. In Figure 2.3 we show a 10-level Haar MRA of the signal f shown in Figure 2.1(a). This signal has 210 values so 10 levels of MRA are possible. On the top of Figure 2.3(a), the graph of A10 is shown; it consists of a single value repeated 210 times. This value is the average of all 210 values of the signal f . The graph directly below it is of the signal A9 which equals A10 plus the details in D10 . Each successive averaged signal is shown, from A10 through A1 . By successively adding on details, the full signal in Figure 2.1(a) is systematically constructed.

2.5

Signal compression

In Section 2.2 we saw that the Haar transform can be used to localize the energy of a signal into a shorter subsignal. In this section we show how this redistribution of energy can be used to compress audio signals. By compressing an audio signal we mean converting the signal data into a new format that requires less bits to transmit. When we use the term, audio signal, we are speaking somewhat loosely. Many of the signals we have in mind are indeed the result of taking discrete samples of a sound signal—as in the data in a computer audio ﬁle, or on a compact disc—but the techniques developed here

0

0.25

0.5

(a)

0.75

1

0

0.25

0.5

0.75

1

(b)

FIGURE 2.3 Haar MRA of the signal in Figure 2.1(a). The graphs are of the ten averaged signals A10 through A1 . Beginning with signal A10 on the top left down to A6 on the bottom left, then A5 on the top right down to A1 on the bottom right.

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

22

also apply to digital data transmissions and to other digital signals, such as digitized electrocardiograms or digitized electroencephalograms. There are two basic categories of compression techniques. The ﬁrst category is lossless compression. Lossless compression methods achieve completely error free decompression of the original signal. Typical lossless methods are Huﬀman compression, LZW compression, arithmetic compression, or run-length compression. These techniques are used in popular lossless compression programs, such as the kind that produce .zip ﬁles. Unfortunately, the compression ratios that can be obtained with lossless methods are rarely more than 2:1 for audio ﬁles consisting of music or speech. The second category is lossy compression. A lossy compression method is one which produces inaccuracies in the decompressed signal. Lossy techniques are used when these inaccuracies are so small as to be imperceptible. The advantage of lossy techniques over lossless ones is that much higher compression ratios can be attained. With wavelet compression methods, if we are willing to accept slight but imperceptible inaccuracies in the decompressed signal, then we can obtain compression ratios of 10:1, or 20:1, or as high as 50:1 or even 100:1. In order to illustrate the general principles of wavelet compression of signals, we shall examine, in a somewhat simpliﬁed way, how the Haar wavelet transform can be used to compress some test signals. For example, Signal 1 in Figure 2.4(a) can be very eﬀectively compressed using the Haar transform. Although Signal 1 is not a very representative audio signal, it is representative of a portion of a digital data transmission. This signal has 1024 values equally spaced over the time interval [0, 20). Most of these values are constant over long stretches, and that is the principal reason that Signal 1 can be compressed eﬀectively with the Haar transform. Signal 2 in Figure 2.5(a), however, will not compress nearly so well; this signal requires the more sophisticated wavelet transforms described in the next chapter. The basic steps for wavelet compression are as follows: Method of Wavelet Transform Compression Step 1. Perform a wavelet transform of the signal. Step 2. Set equal to 0 all values of the wavelet transform which are

insigniﬁcant, i.e., which lie below some threshold value. Step 3. Transmit only the signiﬁcant, non-zero values of the transform obtained from Step 2. This should be a much smaller data set than the original signal. Step 4. At the receiving end, perform the inverse wavelet transform of the data transmitted in Step 3, assigning zero values to the insigniﬁcant values which were not transmitted. This decompression step produces an approximation of the original signal. © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.5 Signal compression

0

5

10

15

4

4

2

2

0

0

−2

−2

−4 20

(a)

0

23

0

5

10

15

−4 20

(b)

5

10

(c)

15

1.01

4

1.005

2

1

0

0.995

−2

0.99 20

0

5

10

15

−4 20

(d)

FIGURE 2.4 (a) Signal 1. (b) 10-level Haar transform of Signal 1. (c) Energy map of Haar transform. (d) 20:1 compression of Signal 1, 100% of energy.

In this chapter we shall illustrate this method using the Haar wavelet transform. This initial discussion will be signiﬁcantly deepened and generalized in the next chapter when we discuss this method of compression in terms of various Daubechies wavelet transforms. Let’s now examine a Haar wavelet transform compression of Signal 1. We begin with Step 1. Since Signal 1 consists of 1024 = 210 values, we can perform 10 levels of the Haar transform. This 10-level Haar transform is shown in Figure 2.4(b). Notice how a large portion of the Haar transform’s values are 0, or very near 0, in magnitude. This fact provides the fundamental basis for performing an eﬀective compression. In order to choose a threshold value for Step 2, we proceed as follows. First, we arrange the magnitudes of the values of the Haar transform so that they are in decreasing order: L1 ≥ L2 ≥ L3 ≥ · · · ≥ LN where L1 is the largest absolute value of the Haar transform, L2 is the next © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

24

2. Haar wavelets

largest, etc. (In the event of a tie, we just leave those magnitudes in their original order.) We then compute the cumulative energy proﬁle of this new signal: 2 L1 L21 + L22 L21 + L22 + L23 , , , ...,1 . Ef Ef Ef For Signal 1, we show a graph of this energy proﬁle—which we refer to as the energy map of the Haar transform—in Figure 2.4(c). Notice that the energy map very quickly reaches its maximum value of 1. In fact, using FAWAV we ﬁnd that L21 + L22 + · · · + L251 = 0.999996. Ef Consequently, if we choose a threshold T that is less than L51 = 0.3536, then the values of the transform that survive this threshold will account for essentially 100% of the energy of Signal 1. We now turn to Step 3. In order to perform Step 3—transmitting only the signiﬁcant transform values—an additional amount of information must be sent which indicates the positions of these signiﬁcant transform values in the thresholded transform. This information is called the signiﬁcance map. The values of this signiﬁcance map are either 1 or 0: a value of 1 if the corresponding transform value survived the thresholding, a value of 0 if it did not. The signiﬁcance map is therefore a string of N bits, where N is the length of the signal. For the case of Signal 1, with a threshold of 0.35, there are only 51 non-zero bits in the signiﬁcance map out of a total of 1024 bits. Therefore, since most of this signiﬁcance map consists of long stretches of zeros, it can be very eﬀectively compressed using one of the lossless compression algorithms mentioned above. This compressed string of bits is then transmitted along with the non-zero values of the thresholded transform. Finally, we arrive at Step 4. At the receiving end, the signiﬁcance map is used to insert zeros in their proper locations in between the non-zero values in the thresholded transform, and then an inverse transform is computed to produce an approximation of the signal. For Signal 1 we show the approximation that results from using a threshold of 0.35 in Figure 2.4(d). This approximation used only 51 transform values; so it represents a compression of Signal 1 by a factor of 1024:51, i.e., a compression factor of 20:1. Since the compressed signal contains nearly 100% of the energy of the original signal, it is a very good approximation. In fact, the maximum error over all values is no more than 1.95 × 10−3 . Life would be simpler if the Haar transform could be used so eﬀectively for all signals. Unfortunately, if we try to use the Haar transform for threshold compression of Signal 2 in Figure 2.5(a), we get poor results. This signal, when played over a computer sound system, produces a sound similar to two low notes played on a clarinet. It has 4096 = 212 values; so we can perform 12 levels of the Haar transform. In Figure 2.5(b) we show a plot of the 12-level Haar transform of Signal 2. It is clear from this plot that a large fraction of © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.5 Signal compression

0

0.5

1

1.5

2

4

1

2

0

0

−1

−2

−2 2

(a)

0

25

0

0.5

1

1.5

−4 2

(b)

0.5

1

(c)

1.5

1.01

2

1.005

1

1

0

0.995

−1

0.99 2

0

0.5

1

1.5

−2 2

(d)

FIGURE 2.5 (a) Signal 2. (b) 12-level Haar transform of Signal 2. (c) Energy map of Haar transform. (d) 10:1 compression of Signal 2, 99.6% of energy.

the Haar transform values have signiﬁcant magnitude, signiﬁcant enough that they are visible in the graph. In fact, the energy map for the transform of Signal 2, shown in Figure 2.5(c), exhibits a much slower increase towards 1 in comparison with the energy map for the transform of Signal 1. Therefore, many more transform values are needed in order to capture a high percentage of the energy of Signal 2. In Figure 2.5(d), we show a 10:1 compression of Signal 2 which captures 99.6% of the energy of Signal 2. Comparing this compression with the original signal we see that it is a fairly poor approximation. Many of the signal values are clumped together in the compressed signal, producing a very ragged or jumpy approximation of the original signal. When this compressed version is played on a computer sound system, it produces a screechy “metallic” version of the two clarinet notes, which is not a very satisfying result. As a rule of thumb, we must capture at least 99.99% of the energy of the signal in order to produce an acceptable approximation, i.e., an approximation that is not perceptually diﬀerent from the original. Achieving this accurate an approximation for Signal 2 requires at least 1782 transform © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

26

values. Because Signal 2 itself has 4096 values, this is a compression ratio of only about 2.3:1, which is not very high. We shall see in the next chapter that Signal 2 can be compressed very eﬀectively, but we shall need more high powered wavelet transforms to do it.

2.5.1

A note on quantization

The most serious oversimpliﬁcation that we made in the discussion above is that we ignored the issue known as quantization. The term quantization is used whenever it is necessary to take into account the ﬁnite precision of numerical data handled by digital methods. For example, the numerical data used to generate the graphs of Signals 1 and 2 above were IEEE double precision numbers that use 8 bytes = 64 bits for each number. In order to compress this data even further, we can represent the wavelet transform coeﬃcients using less bits. We shall address this issue of quantization in the next chapter when we look again at the problem of compression.

2.6

Removing noise

In this section we shall begin our treatment of one of the most important aspects of signal processing, the removal of noise from signals. Our discussion in this section will introduce the fundamental ideas involved in the context of the Haar transform. In the next chapter we shall considerably deepen and generalize these ideas, in the context of the more powerful Daubechies wavelet transforms. When a signal is received after transmission over some distance, it is frequently contaminated by noise. The term noise refers to any undesired change that has altered the values of the original signal. The simplest model for acquisition of noise by a signal is additive noise, which has the form (contaminated signal) = (original signal) + (noise).

(2.30)

We shall represent this equation in a more compact way as f =s+n

(2.31)

where f is the contaminated signal, s is the original signal, and n is the noise. There are several kinds of noise. A few of the commonly encountered types are the following: 1. Random noise. The noise signal is highly oscillatory, its values switching rapidly between values above and below an average, or mean, value. For simplicity, we shall examine random noise with a mean value of 0. 2. Pop noise. This type of noise is heard on old analog recordings obtained from phonograph records. The noise is perceived as randomly occurring, isolated “pops.” As a model for this type of noise we add a few non-zero values to the original signal at isolated locations. © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.6 Removing noise

27

3. Localized random noise. Sometimes the noise appears as in type 1, but only over a short segment or segments of the signal. This occurs when there is a short-lived disturbance in the environment during transmission of the signal. Of course, there can also be noise signals which combine aspects of each of these types. In this section we shall examine only the ﬁrst type of noise, random noise. Other types will be considered later. Our approach will be similar to how we treated compression in the last section; we shall examine how noise removal is performed on two test signals using the Haar transform. For the ﬁrst test signal, the Haar transform is used very eﬀectively for removing the noise. For the second signal, however, the Haar transform performs poorly, and we shall need to use more sophisticated wavelet transforms to remove the noise from this signal. The essential principles, however, underlying these more sophisticated wavelet methods are the same principles we describe here for the Haar transform. Here is our basic method for removing random noise. Threshold Method of Wavelet Denoising Suppose that the contaminated signal f equals the transmitted signal s plus the noise signal n. Also suppose that these two conditions hold: 1. The energy of s is captured, to a high percentage, by transform values whose magnitudes are all greater than a threshold Ts > 0. 2. The noise signal’s transform values all have magnitudes which lie below a noise threshold Tn satisfying Tn < Ts . Then the noise in f can be removed by thresholding its transform: All values of its transform whose magnitudes lie below the noise threshold Tn are set to 0 and an inverse transform is performed, providing a good approximation of f . Let’s see how this method applies to Signal A shown in Figure 2.6(a). This signal was obtained by adding random noise, whose values oscillate between ±0.1 with a mean of zero, to Signal 1 shown in Figure 2.4(a). In this case, Signal 1 is the original signal and Signal A is the contaminated signal. As we saw in the last section, the energy of Signal 1 is captured very eﬀectively by the relatively few transform values whose magnitudes lie above a threshold of 0.35. So we set Ts equal to 0.35, and condition 1 in the Denoising Method is satisﬁed. Now as for condition 2, look at the 10-level Haar transform of Signal A shown in Figure 2.6(b). Comparing this Haar transform with the Haar transform of Signal 1 in Figure 2.4(b), it is clear that the added noise has contributed a large number of small magnitude values to the transform of Signal A, while the high-energy transform values of Signal 1 are plainly visible © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

28

(although slightly altered by the addition of noise). Therefore, we can satisfy condition 2 and eliminate the noise if we choose a noise threshold of, say, Tn = 0.2. This is indicated by the two horizontal lines shown in Figure 2.6(b); all transform values lying between ±0.2 are set to 0, producing the thresholded transform shown in Figure 2.6(c). Comparing Figure 2.6(c) with Figure 2.4(b) we see that the thresholded Haar transform of the contaminated signal is a close match to the Haar transform of the original signal. Consequently, after performing an inverse transform on this thresholded signal, we obtain a denoised signal that is a close match to the original signal. This denoised signal is shown in Figure 2.6(d), and it is clearly a good approximation to Signal 1, especially considering how much noise was originally present in Signal A.

0

5

10

15

4

4

2

2

0

0

−2

−2

−4 20

(a)

0

0

5

10

15

−4 20

(b)

5

10

(c)

15

4

4

2

2

0

0

−2

−2

−4 20

0

5

10

15

−4 20

(d)

FIGURE 2.6 (a) Signal A, 210 values. (b) 10-level Haar transform of Signal A. The two horizontal lines are at ±0.2 where 0.2 is a denoising threshold. (c) Thresholded transform. (d) Denoised signal.

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.6 Removing noise

2.6.1

29

RMS Error

The eﬀectiveness of noise removal can be quantitatively measured in the following way. The Root Mean Square Error (RMS Error) of the contaminated signal f compared with the original signal s is deﬁned to be (f1 − s1 )2 + (f2 − s2 )2 + · · · + (fN − sN )2 (RMS Error) = . (2.32) N For example, for Signal A the RMS Error between it and Signal 1 is 0.057. After denoising, the RMS Error between the denoised signal and Signal 1 is 0.011, a ﬁve-fold reduction in error. This gives quantitative evidence for the eﬀectiveness of the denoising of Signal A. Summarizing this example, we can say that the denoising was eﬀective for two reasons: (1) the transform was able to compress the energy of the original signal into a few high-energy values, and (2) the added noise was transformed into low-energy values. Consequently, the high-energy transform values from the original signal stood out clearly from the low-energy noise transform values which could then be eliminated by thresholding. Unfortunately, denoising with the Haar transform is not always so eﬀective. Consider, for example, Signal B shown in Figure 2.7(a). This signal consists of Signal 2, shown in Figure 2.5(a), with random noise added. We view Signal 2 as the original signal and Signal B as the contaminated signal. As with the ﬁrst case considered above, the random noise has values that oscillate between ±0.1 with a mean of zero. In this case, however, we saw in the last section that it takes a relatively large number of transform values to capture the energy in Signal 2. Most of these transform values are of low energy, and it takes many of them to produce a good approximation of Signal 2. When the random noise is added to Signal 2, then the Haar transform, just like in the previous case, produces many small transform values which lie below a noise threshold. This is illustrated in Figure 2.7(b) where we show the 12-level Haar transform of Signal B. As can be seen by comparing Figure 2.7(b) with Figure 2.5(b), the small transform values that come from the noise obscure most of the small magnitude transform values that result from the original signal. Consequently, when a thresholding is done to remove the noise, as indicated by the horizontal lines in Figure 2.7(b), this removes many of the transform values of the original signal which are needed for an accurate approximation. This can be veriﬁed by comparing the thresholded transform shown in Figure 2.7(c) with the original signal’s transform in Figure 2.5(b). In Figure 2.7(d) we show the denoised signal obtained by inverse transforming the thresholded signal. This denoised signal is clearly an unsatisfactory approximation of the original signal. By computing RMS Errors, we can quantify this judgment. The RMS Error between Signal B and Signal 2 is 0.057, while the RMS Error between the denoised signal and Signal 2 is 0.035. This shows that the error after denoising is almost two-thirds as great as the original error. For this second test case, we can say that the denoising was not eﬀective © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

30

0

0.5

1

1.5

2

2

1

1

0

0

−1

−1

−2 2

(a)

0

0

0.5

1

1.5

−2 2

(b)

0.5

1

1.5

(c)

2

2

1

1

0

0

−1

−1

−2 2

0

0.5

1

1.5

−2 2

(d)

FIGURE 2.7 (a) Signal B, 212 values. (b) 12-level Haar transform of Signal B. The two horizontal lines are at ±0.2 where 0.2 is the denoising threshold. (c) Thresholded transform. (d) Denoised signal.

because the transform could not compress the energy of the original signal into a few high-energy values lying above the noise threshold. We shall see in the next chapter that more sophisticated wavelet transforms can attain the desired compression and thus perform nearly as well at denoising Signal B as the Haar transform did for Signal A.

2.7

Notes and references

More material on the Haar transform and its applications can be found in [3]. Besides the lossy compression method described in this chapter, the Haar transform has also played a role in lossless image compression; see [4] or [5]. 1. R.W. Hamming. (1987). Dover, New York, NY.

Numerical Methods for Scientists and Engineers.

2. L.N. Trefethen. (1998). Maxims about numerical mathematics, computers, science, and life. SIAM News, Vol. 31, No. 1. © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.8 Examples and exercises

31

3. K.R. Rao. (1982). Fast Transforms: Algorithms, Analyses, Applications. Academic Press, New York, NY. 4. M. Rabbini and P.W. Jones. (1991). Digital Image Compression Techniques. SPIE Press, Bellingham, WA. 5. V.K. Heer and H-E Reinfelder. (1990). Comparison of reversible methods for data compression. In Medical Imaging IV, 354–365. Proceedings SPIE, No. 1233.

2.8

Examples and exercises

Note: Computer exercises, designed for FAWAV, are indicated by a superscript c. For example, problem 2.1.5 is a computer exercise. A subscript s means that a solution is provided. For instance, there is a solution provided for problem 2.1.1(a). Solutions are in Appendix B.

Section 2.1 Example 2.1.A

Find the ﬁrst level Haar transform of f = (2, 2, 2, 4, 4, 4).

Solution. The average of the ﬁrst pair of values is 2, the average of the second pair of values is 3, these √ and the average of the √ third √ pair √ of values is 4. Multiplying 1 averages by 2,√we obtain a1 = (2 2, 3 2,√4 2). To√compute √ d , we ﬁnd that 2 = −√ 2, √ and d3 = d1 =√(f1 − f2 )/ 2 = 0, and d2 = (f3 − f4 )/ 2 = −2/ √ √ (f5 − f6 )/ 2 = 0. Thus the ﬁrst level Haar transform of f is (2 2, 3 2, 4 2 | 0, − 2, 0). Example 2.1.B For the signal f = (2, 2, 2, 4, 4, 8), compute an approximate signal e f by inverse transforming the compressed Haar transform (a1 | 0, 0, 0) obtained by setting all the ﬂuctuation values equal to zero. Find the largest error between each value of f and e f. √ √ √ Solution. We ﬁnd that a1 = (2 2, 3 2, 6 2). Therefore, the inverse Haar transform produces e f = (2, 2, 3, 3, 6, 6). The largest error between each value of f and e f is 2. Example 2.1.C [Figure 2.1] To create Figure 2.1 you do the following. First, choose New 1-dim from FAWAV’s menu, and then choose Graph/Plot. Plot the formula3 20 x^2 (1-x)^4 cos(12 pi x) over the interval of type [0, L] with L = 1 using 1024 points. That produces the graph in Figure 2.1(a) (after selecting View/Display style and choosing Blank for the Grid style and Lines for the Plot style). To produce Figure 2.1(b), select Transforms/Wavelet and choose Haar as the Wavelet type with 1 entered for the Levels value. After plotting the transform, change to Lines as the plot style to get the graph shown in the ﬁgure. 2.1.1

Compute the ﬁrst trend and ﬁrst ﬂuctuation for the following signals:

(a)s

f = (2, 4, 6, 6, 4, 2)

(b)

f = (−1, 1, 2, −2, 4, −4, 2, 2)

formula is saved in the archive Primer Formulas.zip which can be downloaded from the book’s website. After extracting the formulas from this archive, you can retrieve them using the Load button under the text box in the graphing procedure. 3 This

© 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2. Haar wavelets

32 (c)s

f = (1, 2, 3, 3, 2, 2, 1, 1)

(d)

f = (2, 2, 4, 4, 6, 6, 8, 8)

2.1.2 Given the following Haar transformed signals, ﬁnd the original signals f that correspond to them using Formula (2.6). √ √ √ √ √ √ (2 2, − 2, 3 2, − 2 | 0, 2, 0, 2) (a)s √ √ √ √ √ √ √ (b) (4 2, 3 2, − 2, 2 2 | 2, − 2, 0, 2 2) √ √ √ √ √ √ (3 2, 2 2, 2 2, 0 | 2 2, − 2, 2, 0) (c)s √ √ √ √ √ √ √ √ (d) (4 2, 5 2, 7 2, −4 2 | 2, 2 2, −2 2, 2) 2.1.3 For each of the signals f given below, compute an approximate signal e f by inverse transforming the compressed Haar transform (a1 | 0, . . . , 0) obtained by setting all the ﬂuctuation values equal to zero. In each case, ﬁnd the largest error between each value of f and e f. (a)s

f = (2, 2, 3, 3, 4, 5, 6, 6)

(b)

f = (1, 2, 3, 3, 2, 1)

(c)s

f = (2, −2, −2, −2, −2, 0, 2, 2)

(d)

f = (4, 4, 4, −1, −1, 1, 2, 2, 4, 6)

2.1.4 Consider again problem 2.1.3 above. When will there be a diﬀerence between a value of f and a value of the approximate signal e f , and when will the two signals’ values be the same? 2.1.5c Plot 1-level Haar transforms of the following functions—sampled uniformly over [0, 1] using 1024 points. (a)

f (x) = x2 (1 − x)

(b)s

f (x) = x4 (1 − x)6 cos(64πx)

(c)

(0.2 < x < 0.3) − 3(0.4 < x < 0.5) + 2(0.5 < x < 0.8)

(d)

f (x) = sgn(sin 12πx)

Section 2.2 Example 2.2.A For the signal f = (2, 2, 4, 6, 8, 10), ﬁnd the energies of its trend and ﬂuctuation subsignals and show that their sum equals the energy of f . √ √ √ √ √ Solution. The trend is a1 = (2 2, 5 2, 9 2) and ﬂuctuation is d1 = (0, − 2, − 2). The trend energy Ea1 is 8 + 50 + 162 = 220, and the ﬂuctuation energy is Ed1 = 0 + 2 + 2 = 4. Their sum is 224 and the energy of f is 4 + 4 + 16 + 36 + 64 + 100 = 224 so they are equal. Example 2.2.B Compute the percentage of compaction of the energy of f = (2, 2, 4, 6, 8, 10) by the 1-level Haar transform, in terms of Ea1 /Ef . Solution. In Example 2.2.A, we found that Ea1 = 220 and that Ef = 224. Therefore, Ea1 /Ef = 220/224 = 0.982 . . . © 2008 by Taylor & Francis Group, LLC

i

i i

i

i

i

i

i

2.8 Examples and exercises

33

Example 2.2.C [Figure 2.2] To graph Figure 2.2(a), you plot (after selecting Edit/Points used and selecting 4096 as the number of points): 50x^2(1-x)^6cos(12pi x)(0<x

Our partners will collect data and use cookies for ad personalization and measurement. Learn how we and our ad partner Google, collect and use data. Agree & close