An Introduction To Digital Signal Processing: Edmund M-K. Lai

An Introduction to Digital Signal Processing
Edmund M-K. Lai 17th November 2003
Contents
Preface 1 Introduction 1.1 1.2 1.3 Advantages of Processing Signals Digitally . . . . . . . . . . . . . . . . Denition of Some Terms . . . . . . . . . . . . . . . . . . . . . . . . . . DSP Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 1.3.2 1.3.3 1.4 1.4.1 Sample-by-sample Processing . . . . . . . . . . . . . . . . . . . Block Processing . . . . . . . . . . . . . . . . . . . . . . . . . . Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Speech and Audio Processing . . . . . . . . . . . . . . . . . . . 1.4.1.1 1.4.1.2 1.4.1.3 1.4.1.4 1.4.2 1.4.2.1 1.4.2.2 1.4.2.3 1.4.3 Speech Coding . . . . . . . . . . . . . . . . . . . . . . Linear Prediction . . . . . . . . . . . . . . . . . . . . Speech Synthesis . . . . . . . . . . . . . . . . . . . . Speech Recognition . . . . . . . . . . . . . . . . . . . Image Enhancement . . . . . . . . . . . . . . . . . . . xi 1 2 2 4 4 4 5 5 5 5 7 9 9 9 9
Application Areas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Image and Video Processing . . . . . . . . . . . . . . . . . . . .
Image Restoration . . . . . . . . . . . . . . . . . . . . 10 Image Compression and Coding . . . . . . . . . . . . 10
Adaptive Filtering . . . . . . . . . . . . . . . . . . . . . . . . . 11
iv 1.4.3.1 1.4.3.2 1.4.3.3 1.4.4 1.4.5 1.4.6 1.5 2
CONTENTS Noise Cancellation . . . . . . . . . . . . . . . . . . . . 11 Echo Cancellation . . . . . . . . . . . . . . . . . . . . 11 Channel Equalization . . . . . . . . . . . . . . . . . . 12
Control Applications . . . . . . . . . . . . . . . . . . . . . . . . 13 Sensor or Antenna Array Processing . . . . . . . . . . . . . . . 15
Digital Communication Transceivers . . . . . . . . . . . . . . . 16
Brief Overview of Chapters . . . . . . . . . . . . . . . . . . . . . . . . . 16 19
Discrete-time Signals and Systems 2.1
Discrete-time signals . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.1.1 Some Elementary Sequences . . . . . . . . . . . . . . . . . . . . 20 2.1.1.1 2.1.1.2 2.1.1.3 2.1.1.4 2.1.1.5 2.1.2 2.1.3 Unit Impulse Sequence . . . . . . . . . . . . . . . . . 20 Unit Step Sequence . . . . . . . . . . . . . . . . . . . 20 Sinusoidal Sequences . . . . . . . . . . . . . . . . . . 21 Complex Exponential Sequences . . . . . . . . . . . . 21 Random Sequences . . . . . . . . . . . . . . . . . . . 22
Types of Sequences . . . . . . . . . . . . . . . . . . . . . . . . . 22 Some Basic Operations on Sequences . . . . . . . . . . . . . . . 27
2.2
Discrete-time Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 2.2.1 2.2.2 Classication of Systems . . . . . . . . . . . . . . . . . . . . . . 30 Linear Shift-Invariant Systems . . . . . . . . . . . . . . . . . . . 34 2.2.2.1 2.2.2.2 2.2.2.3 2.2.2.4 2.2.3 Linear Convolution . . . . . . . . . . . . . . . . . . . 34 Properties of Linear Convolution . . . . . . . . . . . . 38 Condition for Stability . . . . . . . . . . . . . . . . . . 38 Condition for Causality . . . . . . . . . . . . . . . . . 39
FIR and IIR Systems . . . . . . . . . . . . . . . . . . . . . . . . 40
2.3
Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
CONTENTS 3 Discrete-time Fourier Analysis 3.1 3.2
v 43
Continuous-time Fourier Analysis . . . . . . . . . . . . . . . . . . . . . 43 Discrete-time Fourier Series . . . . . . . . . . . . . . . . . . . . . . . . 44 3.2.1 3.2.2 3.2.3 Average Power . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 Normalized Frequency . . . . . . . . . . . . . . . . . . . . . . . 47 Power Density Spectrum . . . . . . . . . . . . . . . . . . . . . . 48
3.3
Discrete-time Fourier Transform . . . . . . . . . . . . . . . . . . . . . . 49 3.3.1 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.3.1.1 3.3.1.2 3.3.2 3.3.3 Uniform Convergence . . . . . . . . . . . . . . . . . . 50 Mean Square Convergence . . . . . . . . . . . . . . . 50
Energy Density Spectrum . . . . . . . . . . . . . . . . . . . . . 52 Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.3.3.1 3.3.3.2 3.3.3.3 3.3.3.4 3.3.3.5 3.3.3.6 3.3.3.7 Periodicity . . . . . . . . . . . . . . . . . . . . . . . . 53 Symmetry . . . . . . . . . . . . . . . . . . . . . . . . 53 Linearity . . . . . . . . . . . . . . . . . . . . . . . . . 54 Time Shifting . . . . . . . . . . . . . . . . . . . . . . 54 Time Reversal . . . . . . . . . . . . . . . . . . . . . . 55 Frequency Shifting . . . . . . . . . . . . . . . . . . . 55 Linear Convolution . . . . . . . . . . . . . . . . . . . 55
3.3.4 3.3.5 3.4
Inverse DTFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 Autocorrelation Sequence . . . . . . . . . . . . . . . . . . . . . 56
Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 59
4 Discrete Fourier Transform and Its Computation 4.1
Discrete Fourier Transform . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.1.1 4.1.2 Sampling the DTFT . . . . . . . . . . . . . . . . . . . . . . . . 59 Inverse DFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
vi 4.1.3 4.1.4 4.1.5 4.1.6 4.1.7 4.1.8 4.1.9 4.2
CONTENTS Matrix Form . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 Zero-padding . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 Circular Shift . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 Circular Convolution . . . . . . . . . . . . . . . . . . . . . . . . 67 Linear Convolution Using DFT . . . . . . . . . . . . . . . . . . 69 Computational Complexity . . . . . . . . . . . . . . . . . . . . . 75
The Fast Fourier Transform . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.2.1 4.2.2 4.2.3 Radix-2 Decimation-in-time Algorithm . . . . . . . . . . . . . . 77 Radix-2 Decimation-in-frequency Algorithm . . . . . . . . . . . 81 Other Fast Algorithms . . . . . . . . . . . . . . . . . . . . . . . 84
4.3
Spectral Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.3.1 4.3.2 4.3.3 4.3.4 4.3.5 Spectral Leakage . . . . . . . . . . . . . . . . . . . . . . . . . . 84 Rectangular Windowing . . . . . . . . . . . . . . . . . . . . . . 86 Frequency Resolution . . . . . . . . . . . . . . . . . . . . . . . 89 Computational Frequency Resolution . . . . . . . . . . . . . . . 92 Non-stationary Signals . . . . . . . . . . . . . . . . . . . . . . . 93
4.4 5
Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 97
Sampling, Reconstruction and Quantization 5.1
The Sampling Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 5.1.1 5.1.2 5.1.3 5.1.4 Spectra of Sampled Signals . . . . . . . . . . . . . . . . . . . . 99 Aliasing and The Uniform Sampling Theorem . . . . . . . . . . 100 Anti-aliasing Prelters . . . . . . . . . . . . . . . . . . . . . . . 107 Practical Limits on Sampling Rates . . . . . . . . . . . . . . . . 108
5.2
Reconstruction of Analog Signals . . . . . . . . . . . . . . . . . . . . . 108 5.2.1 Ideal Reconstructor . . . . . . . . . . . . . . . . . . . . . . . . . 109
CONTENTS 5.2.2 5.2.3 5.3
vii Staircase Reconstructor . . . . . . . . . . . . . . . . . . . . . . . 111 Anti-image Postlters . . . . . . . . . . . . . . . . . . . . . . . 111
Quantization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 5.3.1 5.3.2 5.3.3 5.3.4 Sample-and-Hold . . . . . . . . . . . . . . . . . . . . . . . . . . 114 Uniform Quantization . . . . . . . . . . . . . . . . . . . . . . . 115 Quantization Noise . . . . . . . . . . . . . . . . . . . . . . . . . 116 Signal-to-Quantization Noise Ratio . . . . . . . . . . . . . . . . 117
5.4
Sampling Rate and Quantization Tradeoff . . . . . . . . . . . . . . . . . 118 5.4.1 5.4.2 5.4.3 5.4.4 Power Spectral Density . . . . . . . . . . . . . . . . . . . . . . . 118 Oversampling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 Noise Shaping . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 Non-uniform Quantization . . . . . . . . . . . . . . . . . . . . . 122 5.4.4.1 Companding . . . . . . . . . . . . . . . . . . . . . . . 123
5.5 5.6 6
Dithering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 131
Digital Filters and Their Realizations 6.1
A Review of Z Transform . . . . . . . . . . . . . . . . . . . . . . . . . . 131 6.1.1 6.1.2 6.1.3 Denition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131 Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 Inverse Z Transform . . . . . . . . . . . . . . . . . . . . . . . . 133
6.2
System Transfer Functions . . . . . . . . . . . . . . . . . . . . . . . . . 134 6.2.1 6.2.2 Poles and Zeros . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 Frequency Response . . . . . . . . . . . . . . . . . . . . . . . . 134
6.3
Digital Filter Realization . . . . . . . . . . . . . . . . . . . . . . . . . . 137 6.3.1 6.3.2 Direct Form I . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138 Canonical Form or Direct Form II . . . . . . . . . . . . . . . . . 139
viii 6.3.3 6.3.4 6.3.5 6.3.6 6.3.7 6.4 6.4.1 6.4.2 6.4.3 6.4.4 6.4.5 6.4.6 6.5 7
CONTENTS Transposed Structures . . . . . . . . . . . . . . . . . . . . . . . 143 FIR Filters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 Cascade Realization . . . . . . . . . . . . . . . . . . . . . . . . 146 Parallel Realization . . . . . . . . . . . . . . . . . . . . . . . . . 149 Circular Buffers . . . . . . . . . . . . . . . . . . . . . . . . . . 150 Fixed-point Representation . . . . . . . . . . . . . . . . . . . . . 155 Floating-point Representation . . . . . . . . . . . . . . . . . . . 157 Coefcient Quantization . . . . . . . . . . . . . . . . . . . . . . 158 Round-off Noise . . . . . . . . . . . . . . . . . . . . . . . . . . 162 Limit Cycles . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165 Arithmetic Overow . . . . . . . . . . . . . . . . . . . . . . . . 167
Finite Wordlength Effects . . . . . . . . . . . . . . . . . . . . . . . . . . 154
Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169 171
FIR Filter Design 7.1 7.2
Design Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171 Window Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176 7.2.1 7.2.2 Rectangular Window . . . . . . . . . . . . . . . . . . . . . . . . 177 Kaiser Window . . . . . . . . . . . . . . . . . . . . . . . . . . . 180
7.3 7.4 7.5 7.6 8
Highpass Filters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181 Bandpass Filters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 182 Frequency Sampling Method . . . . . . . . . . . . . . . . . . . . . . . . 183 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184 187
IIR Filter Design 8.1 8.2
Design Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 188 Some Classic Analog Filters . . . . . . . . . . . . . . . . . . . . . . . . 188 8.2.1 Butterworth Filters . . . . . . . . . . . . . . . . . . . . . . . . . 188
CONTENTS 8.2.2 8.3
ix Chebyshev Filters . . . . . . . . . . . . . . . . . . . . . . . . . . 192
Bilinear Transformation . . . . . . . . . . . . . . . . . . . . . . . . . . . 196 8.3.1 8.3.2 Frequency Warping . . . . . . . . . . . . . . . . . . . . . . . . . 196 Frequency Transformation . . . . . . . . . . . . . . . . . . . . . 199
8.4 8.5
Impulse Invariant Method . . . . . . . . . . . . . . . . . . . . . . . . . . 202 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 206 207 209
Bibliography Index
CONTENTS
Preface
Many excellent textbooks have been written on the subject of Digital Signal Processing (DSP). See the bibliography at the end of the book for examples. They typically cover a wide range of topics and in much detail. Since DSP is such a broad area, there are very few DSP textbooks that are less than 500 pages in length. Students often nd themselves overwhelmed by the depth and breadth covered by these books. Moreover, the amount of mathematics used gave students the impression that DSP is a very difcult subject. I have been teaching DSP to undergraduates and practising engineers and technicians over a number of years. Often the students both in universities and in the industry do not realize that it is some of the fundamental concepts in DSP that will take time to digest. Once they grasped some of these concepts, the mathematics becomes obvious and simple. The emphasis of this book is on the fundamental concepts and ideas that are important and common to most DSP applications. I have consciously keep the length of the book to a minimum by carefully selecting the topics to be covered. That means of course that there are a number of interesting topics in DSP that are not even touched upon. Readers are encouraged to consult the more comprehensive DSP texts if they are interested in extending their knowledge to these other areas. The fundamentals that they have learnt from this book should enable them to pursue the subject further with emphases based on their areas of interest. The materials in this book are based on the course I give to senior Computer Engineering undergraduates in Nanyang Technological University, Singapore. They are suitable for use in a single semester course for senior undergraduate or rst year graduate students. Readers are assumed to have some basic knowledge of continuous-time signals and systems, and Z transform. A typical mathematics or signal and systems course in the early part of an electrical and computer engineering curriculum will sufce.
A The book is typeset entirely using LTEX and Lyx on a Linux platform. Most of the diagrams are drawn using the Xg program. All these are open source programs that
xii
Preface
are the efforts of thousands of volunteers. I wish to thank them for producing such high quality software for the common good. The remaining graphs and plots are produced using MATLAB R . Many former students of The University of Western Australia, The Chinese University of Hong Kong, Edith Cowan University in Perth, Australia, and Nanyang Technological University, Singapore since the late 1980s have been taught using part or all of the materials in this book. To see them learn and develop an interest in DSP has always been a joy. They in turn inspired me by asking insightful questions and by challenging my presentations. I wish to thank all of them. I also wish to thank my wife, Susie, and my daughters, Rhoda and Sarah, for bring so much joy and purpose to my life. Without them it would not have been worth the effort at all. This book is dedicated to them. Singapore
Chapter 1 Introduction
Digital signal processing (DSP) is concerned with the use of digital processors to analyze, modify and extract information from a signal. Signals are normally physical quantities that can be measured. DSP is a eld which is primarily technology driven. It started from around the mid-1960s when digital computers and digital circuitry became fast enough to process large amount of data efciently. When the term digital is used, often it loosely refers to a nite set of distinct values. This is in contrast to analog which refers to a continuous range of values. In digital signal processing we are concerned with the processing of signals which are discrete in time (sampled) and, in most cases, discrete in amplitude (quantized) as well. In other words, we are primarily dealing with data sequences sequences of numbers. Such discrete (or digital) signals may arise in one of the following two distinct circumstances: 1. The signal may be inherently discrete in time (and/or amplitude); or 2. The signal may be a sampled version of a continuous-time signal. Examples of the rst type of data sequences include monthly sales gures, daily highest/lowest temperatures, stock market indices and students examination marks. Business people, meteorologists, economists, and teachers process these types of data sequences to determine cyclic patterns, trends, and averages. The processing usually involves ltering to remove as much noise as possible so that the pattern of interest will be enhanced or highlighted.
Introduction
Examples of the second type of discrete-time signals can readily be found in many engineering applications. For instance, speech and audio signals are sampled and then encoded for storage or transmission. A compact disc player reads the encoded digital audio signals and reconstructs the continuous-time signals for playback.
1.1
Advantages of Processing Signals Digitally
A typical question one may ask is: why process signals digitally? For the rst type of signals discussed previously, the reason is obvious. If the signals are inherently discrete in time, the most natural way to process them is using digital methods. But for continuoustime signals, we have a choice. Analog signals are processed by analog electronics while digital signals can be processed by computers or microprocessors. Analog methods are potentially faster since the analog circuits process signals as they arrive in real-time, provided the settling time is fast enough. On the other hand, digital techniques are algorithmic in nature. If the computer is fast and the algorithms are efcient, then digital processing can be performed in real-time provided the data rate is slow enough. However, with the speed of digital logic increasing exponentially, the upper limit in data rate that can still be considered as real-time processing is becoming higher and higher. The major advantage of digital signal processing is consistency. For the same signal, the output of a digital process will always be the same. It is not sensitive to offsets and drifts in electronic components. The second main advantage of DSP is that very complex digital logic circuits can be packed onto a single chip, thus reducing the component count and the size and reliability of the system.
1.2
Denition of Some Terms
DSP has its origin in electrical/electronic engineering (EE). Therefore the terminologies used in DSP are typically that of EE. We shall now attempt to explain a few terms that we shall be using throughout the book. Signals we have already started using this term in the previous section. A signal is simply a quantity that we can measure over a period of time. This quantity usually
Denition of Some Terms
changes with time and that is what makes it interesting. Such quantities could be voltage or current. They could also be the pressure, uid level and temperature. Other quantities of interest include nancial indices such as the stock market index. You will be surprised how much of the concepts in DSP has been used to analyse the nancial market. Frequency some signals change slowly over time and others change rapidly. For instance, the (AC) voltage available at our household electrical mains goes up and down like a sine function and they complete one cycle in 50 times or 60 times a second. This signal is said to have a frequency of 50 or 60 Hertz (Hz). Spectrum while some signals consist of only a single frequency, others have a combination of a range of frequencies. If you play a string on the violin, there is a fundamental tone (frequency) corresponding to the musical note that is played. But there are other harmonics (integer multiples of the fundamental frequency) present. This musical sound signal is said to have a spectrum of frequencies. The spectrum is a frequency (domain) representation of the time (domain) signal. The two representations are equivalent. Low-pass Filter lters let a certain range of frequency components of a signal through while rejecting the other frequency components. A low-pass lter lets the "lowfrequency" components through. Low-pass lters have a cutoff frequency below which the frequency components can pass through the lter. For instance, if a signal has two frequency components, say, 10 Hz and 20 Hz. Applying a low-pass lter to this signal with a cutoff frequency of 15 Hz will result in an output signal which has only one frequency component at 10 Hz; the 20 Hz component has been rejected by the lter. Bandpass Filter bandpass lters are similar to low-pass lters in that only a range of frequency components can pass through it intact. This range (the passband) is usually above the DC (zero frequency) and somewhere in the mid-range. For instance, we can have a bandpass lter with a passband between 15 and 25 Hz. Applying this lter to the signal discussed above will result in a signal having only a 20 Hz component. High-pass Filter these lters allows frequency components above a certain frequency (cutoff) to pass through intact, rejecting the ones lower than the cutoff frequency. This should be enough to get us going. New terms will arise from time to time and they will be explained as we come across them.
Introduction
1.3
DSP Systems
DSP systems are discrete-time systems which means that they accept digital signals as input and output digital signals (or information extracted). Digital signals are simply sequences of numbers. The output sequence of samples y[n] is computed from the input sequence of samples x[n] according to some rules which the discrete-time system H denes. There are two main methods by which the output sequence is computed from the input sequence. They are called sample-by-sample processing and block processing respectively. We shall encounter both types of processing in later chapters. Most systems can be implemented with either processing method. The output obtained in both cases should be equivalent if the input and the system are the same.
1.3.1 Sample-by-sample Processing

With the sample-by-sample processing method, normally one output sample is obtained when one input sample is presented to the system. For instance, if the sequence {y0 , y1 , y2 , . . . , yn , . . .} is obtained when the input sequence {x0 , x1 , x2 , . . . , xn , . . .} is presented to the system. The sample y0 appears at the output when the input x0 is available at the input. The sample y1 appears at the output when the input x1 is available at the input, etc. The delay between the input and output for sample-by-sample processing is at most one sample. The processing has to be completed before the next sample appears at the input.
1.3.2
Block Processing
With block processing methods, a block of signal samples are being processing at a time. A block of samples are usually treated as a vector which is transformed to an output vector of samples by the system transformation H. x0 y0 x1 y1 x= x y = y H (1.1) 2 2 . . . . . .
Application Areas
The delay between input and output in this case is dependent on the number of samples in each block. For example, if we use 8 samples per block, then the rst 8 input samples have to be buffered (or collected) before processing can proceed. So the block of 8 output samples will appear at least 8 samples after the rst sample x0 appears. The block computation (according to H) has to be completed before the next block of 8 samples are collected.
1.3.3 Remarks
Both processing methods are extensively used in real applications. We shall encounter DSP algorithms and implementation that uses one or the other. The reader might nd it useful in understanding the algorithms or techniques being discussed by being conscious of which processing method is being used.
1.4
Application Areas
Digital signal processing is being applied to a large range of applications. No attempt is made to include all areas of application here. In fact, new applications are constantly appearing. In this section, we shall try to describe a sufciently broad range of applications so that the reader can get a feel of what DSP is about.
1.4.1 Speech and Audio Processing

An area where DSP has found a lot of application is in speech processing. It is also one of the earliest applications of DSP. Digital speech processing includes three main sub-areas: encoding, synthesis, and recognition.
1.4.1.1
Speech Coding
There is a considerable amount of redundancy in the speech signal. The encoding process removes as much redundancy as possible while retaining an acceptable quality of the remaining signal. Speech coding can be further divided into two areas:
Introduction 1. compression a compact representation of the speech waveform without regard to its meaning. 2. parametrization a model that characterizes the speech in some linguistically or acoustically meaningful form.
The minimum channel bandwidth required for the transmission of an acceptable quality of speech is around 3kHz, with a dynamic range of 72dB. This is normally referred to as telephone quality. Converting into digital form, a sampling rate of 8k samples per second with a 12-bit quantization (212 amplitude levels) is commonly used, resulting in 96k bits per second of data. This data rate can be signicantly reduced without affecting the quality of the reconstructed speech as far as the listener is concerned. We shall briey describe three of them: Companding or Non-uniform Quantization The dynamic range of speech signals is very large. This is due to the fact that voiced sounds such as vowels contains a lot of energy and exhibits wide uctuations in amplitude while unvoiced sounds like fricatives generally have much lower amplitudes. A compander (compressor-expander) compresses the amplitude of the signal at the transmitter end and expands it at the receiver end. The process is illustrated schematically in Figure 1.1. The compressor compresses the large amplitude samples and expands the small amplitude ones while the expander does the opposite. The -law compander (with = 255) is a North American standard. A-law companding with A = 87.56 is a European (CCITT) standard. The difference in performance is minimal. A-law companding gives slightly better performance at high signal levels while -law is better at low levels. Adaptive Differential Quantization At any adequate sampling rate, consecutive samples of the speech signal are highly correlated, except for those sounds that contain a signicant amount of wideband noise. The data rate can be greatly reduced by quantizing the difference between two samples instead. Since the dynamic range will be much reduced by differencing, the number of levels required for the quantizer will also be reduced. The concept of differential quantization can be extended further. Suppose we have an estimate of the value of the current sample based on information from the previous samples, then we can quantize the difference between the current sample and its estimate. If the prediction is accurate enough, this difference will be quite small.
Application Areas
Transmitter Signal x(t) Receiver Recovered Signal x(t)
Compressor
y(t)
y(t)
Expander
y(t) 1
x(t) 1
x(t)
y(t)
Figure 1.1: The companding process Figure 1.2 shows the block diagram of an adaptive differential pulse code modulator (ADPCM). It takes a 64 kbits per second pulse code modulated (PCM) signal and encodes it into 32kbit per second adaptive differential pulse code modulated (ADPCM) signal. 1.4.1.2 Linear Prediction The linear predictive coding method of speech coding is based on a (simplied) model of speech production shown in Figure 1.3. The time-varying digital lter models the vocal tract and is driven by an excitation signal. For voiced speech, this excitation signal is typically a train of scaled unit impulses at pitch frequency. For unvoiced sounds it is random noise. The analysis system (or encoder) estimates the lter coefcients, detects whether the speech is voiced or unvoiced and estimates the pitch frequency if necessary. This is performed for each overlapping section of speech usually around 10 milliseconds in duration. This information is then encoded and transmitted. The receiver reconstructs the speech signal using these parameters based on the speech production model. It is interesting to note that the reconstructed speech is similar to the original perceptually but the physical appearance of the signal is very different. This is an illustration of the redundancies inherent in speech signals.
Introduction
Step Size Adjustment
[n]
x[n]
error signal
Quantizer
Encoder
c[n]
Predictor
Figure 1.2: Block diagram of an adaptive differential pulse code modulator
Impulse Train For voiced speech
Vocal Tract Model Timevarying Filter Gain Synthesized Speech
Random Noise For unvoiced speech
Voiced/Unvoiced Decision
Figure 1.3: A model of speech production
Application Areas 1.4.1.3 Speech Synthesis
The synthesis or generation of speech can be done through the speech production model mentioned above. Although the duplication of the acoustics of the vocal tract can be carried out quite accurately, the excitation model turns out to be more problematic. For synthetic speech to sound natural, it is essential that the correct allophone be produced. Despite the fact that different allophones are perceived as the same sound, if the wrong allophone is selected, the synthesized speech will not sound natural. Translation from phonemes to allophones is usually controlled by a set of rules. The control of timing of a word is also very important. But these rules are beyond the realm of DSP. 1.4.1.4 Speech Recognition
One of the major goal of speech recognition is to provide an alternative interface between human user and machine. Speech recognition systems can either be speaker dependent or independent, and they can either accept isolated utterances or continuous speech. Each system is capable of handling a certain vocabulary. The basic approach to speech recognition is to extract features of the speech signals in the training phase. In the recognition phase, the features extracted from the incoming signal are compared to those that have been stored. Owing to the fact that our voices change with time and the rate at which we speak also varies, speech recognition is a very tough problem. However, there are now commercially available some relatively simple small vocabulary, isolated utterance recognition systems. This comes about after 30 years of research and the advances made in DSP hardware and software.
1.4.2 Image and Video Processing

Image processing involves the processing of signals which are two dimensional. A digital image consists of a two dimensional array of pixel values instead of a one dimensional one for, say, speech signals. We shall briey describe three areas of image processing. 1.4.2.1 Image Enhancement
Image enhancement is used when we need to focus or pick out some important features of an image. For example, we may want to sharpen the image to bring out details such
10
Introduction
as a car licence plate number or some areas of an X-ray lm. In aerial photographs, the edges or lines may need to be enhanced in order to pick out buildings or other objects. Certain spectral components of an image may need to be enhanced in images obtained from telescopes or space probes. In some cases, the contrast may need to be enhanced. While linear ltering may be all that is required for certain types of enhancement, most useful enhancement operations are nonlinear in nature.
1.4.2.2
Image Restoration
Image restoration deals with techniques for reconstructing an image that may have been blurred by sensor or camera motion, and in which additive noise may be present. The blurring process is usually modeled as a linear ltering operation and the problem of image restoration then becomes one of identifying the type of blur and estimating the parameters of the model. The image is then ltered by the inverse of the lter.
1.4.2.3
Image Compression and Coding
The amount of data in a visual image is very large. A simple black-and-white still picture digitized to a 512 512 array of pixels using 8 bits per pixel involves more than 2 million bits of information. In the case of sequences of images such as in video or television images, the amount of data involved will be even greater. Image compression, like speech compression, seeks to reduce the number of bits required to store or transmit the image with either no loss or an acceptable level of loss or distortion. A number of different techniques have been proposed, including prediction or coding in the (spatial) frequency domain. The most successful techniques typically combine several basic methods. Very sophisticated methods have been developed for digital cameras and digital video discs (DVD). Standards have been developed for the coding of both image and video signals for different kinds of applications. For still images, the most common one is JPEG. For high quality motion video, there is MPEG and MPEG-2. MPEG-2 was developed with highdenition television in mind. It is now used in satellite transmission of broadcast quality video signals.
Application Areas
11
1.4.3 Adaptive Filtering

A major advantage of digital processing is its ability of adapt to changing environments. Even though adaptive signal processing is a more advanced topic which we will not cover in this course, we shall describe the basic ideas involved in adaptive signal processing and some of its applications. A basic component in an adaptive digital signal processing system is a digital lter with adjustable lter coefcients a time-varying digital lter. Changing the characteristics of a lter by a change in the coefcient values is a very simple operation in DSP. The adaptation occurs through an algorithm which takes a reference (or desired) signal and an error signal produced by the difference between the current output of the lter and the current input signal. The algorithm adjusts the lter coefcients so that the averaged error is minimized. 1.4.3.1 Noise Cancellation
One example of noise cancellation is the suppression of the maternal ECG component in fetal ECG. The fetal heart rate signal can be obtained from a sensor placed in the abdominal region of the mother. However, this signal is very noisy due to the mothers heartbeat and fetal motion. The idea behind noise cancellation in this case is to take a direct recording of the mothers heartbeat and after ltering of this signal, subtract it off the fetal heart rate signal to get a relatively noise-freee heart rate signal. A schematic diagram of the system is shown in Figure 1.4. There are two inputs: a primary and a reference. The primary signal is of interest but has a noisy interference component which is correlated with the reference signal. The adaptive lter is used to produce an estimate of this interference or noise component which is then subtracted off the primary signal. The lter should be chosen to ensure that the error signal and the reference signal are uncorrelated. 1.4.3.2 Echo Cancellation
Echoes are signals that are identical to the original signals but are attenuated and delayed in time. They are typically generated in long distance telephone communication due to impedance mismatch. Such a mismatch usually occurs at the junction or hybrid between
12
Introduction
Reference signal from chest sensor
Noise
Primary signal from abdominal sensor
Adjustable Filter Filter coefficient adjustment Adaptive Algorithm
Figure 1.4: An adaptive noise cancellation system the local subscriber loop and the long distance loop. As a result of the mismatch, incident electromagnetics waves are reected which sound like echoes to the telephone user. The idea behind echo cancellation is to predict the echo signal values and thus subtract it out. The basic mechanism is illustrated in Figure 1.5. Since the speech signal is constantly changing, the system has to be adaptive. 1.4.3.3 Channel Equalization
Consider the transmission of a signal over a communication channel (e.g. coaxial cable, optical bre, wireless). The signal will be subject to channel noise and dispersion caused, for example, by reection from objects such as buildings in the transmission path. This distorted signals will have to be reconstructed by the receiver. One way to restore the original signal is to pass the received signal through an equalizing lter to undo the dispersion effects. The equalizer should ideally be the inverse of the channel characteristics. However, channel characteristics typically drift in time and so the equalizer (a digital lter) coefcients will need to be adjusted continuously. If the transmission medium is a cable, the drift will occur very slowly. But for wireless channels in mobile communications the channel characteristics change rapidly and the equalizer lter will have to adapt very quickly.
Application Areas
Nearend signal
13
Farend signal
4 wire to 2 wire junction hybrid
Nearend signal and echo of farend signal
Adjustable Filter Filter coefficient adjustment Adaptive Algorithm
Cleaned nearend signal
To farend
Figure 1.5: An adaptive echo cancellation system In order to learn the channel characteristics, the adaptive equalizer operates in a training mode where a pre-determined training signal is transmitted to the receiver. Normal signal transmission has to be regularly interrupted by a brief training session so that the equalizer lter coefcients can be adjusted. Figure 1.6 shows an adaptive equalizer in training mode.
1.4.4
Control Applications
A digital controller is a system used for controlling closed-loop feedback systems as shown in Figure 1.7. The controller implements algebraic algorithms such as lters and compensators to regulate, correct, or change the behaviour of the controlled system. Digital control has the advantage that complex control algorithms are implemented in software rather than specialized hardware. Thus the controller design and its parameters can easily be altered. Furthermore, increased noise immunity is guaranteed and parame-
14
Introduction
Training Signal
Delay
Channel
Adjustable Filter (Equalizer)
Noise
Adaptive Algorithm
Error Signal
Figure 1.6: An adaptive equaliser in training mode
Reference Signal
error signal
Controller
System
Output
Adaptive Algorithm
Parameter Adjustment
Figure 1.7: A digital closed loop control system
Application Areas
15
ter drift is eliminated. Consequently, they tend to be more reliable and at the same time, features reduced size, power, weight and cost. Digital signal processors are very useful for implementing digital controllers since they are typically optimized for digital ltering operations with single instruction arithmetic operations. Furthermore, if the system being controlled changes with time, adaptive control algorithms, similar to adaptive ltering discussed above, can be implemented.
1.4.5 Sensor or Antenna Array Processing

In some applications, a number of spatially distributed sensors are used for receiving signals from some sources. The problem of coherently summing the outputs from these sensors is known as beamforming. Beyond the directivity provided by an individual sensor, a beamformer permits one to listen preferentially to wavefronts propagating from one direction over another. Thus a beamformer implements a spatial lter. Applications of beamforming can be found in seismology, underwater acoustics, biomedical engineering, radio communication systems and astronomy. In cellular mobile communication systems, smart antennas (an antenna array with digitally steerable beams) are being used to increase user capacity and expand geographic coverage. In order to increase capacity, an array which can increase the carrier to interference ratio (C/I) at both the base station and the mobile terminal is required. There are three approaches to maximizing C/I with an antenna array. 1. The rst one is to create higher gain on the antenna in the intended direction using antenna aperture. This is done by combining the outputs of each individual antenna to create aperture. 2. The second approach is the mitigation of multipath fading. In mobile communication, fast fading induced by multipath propagation requires an additional link margin of 8dB. This margin can be recovered by removing the destructive multipath effects. 3. The third approach is the identication and nulling of interferers. It is not difcult for a digital beamformer to create sharp nulls, removing the effects of interference. Direction of arrival estimation can also be performed using sensor arrays. In the simplest conguration, signals are received at two spatially separated sensors with one signal being an attenuated, delayed and noisy version of the other. If the distance between the
16
Introduction
sensors is known, and the signal velocity is known, then the direction of arrival can be estimated. If the direction does not change, or changes very slowly with time, then it can be determined by cross-correlating the two signals and nding the global maximum of the cross-correlation function. If the direction changes rapidly, then an adaptive algorithm is needed.
1.4.6 Digital Communication Transceivers

One of the most exciting applications of DSP is in the design and implementation of digital communication equipment. Throughout the 1970s and 80s radio systems migrated from analog to digital in almost every aspect, from system control to source and channel coding to hardware technology. A new architecture known generally as software radio is emerging. This technology liberates radio-based services from dependence on hardwired characteristics such as frequency band, channel bandwidth, and channel coding. The software radio architecture centres on the use of wideband analog-to-digital and digital-to-analog converters that are placed as close to the antenna as possible. Since the signal is being digitized earlier in the system, as much radio functionality as possible can be dened and implemented in software. Thus the hardware is relatively simple and functions are software dened as illustrated in Figure 1.8. Software denable channel modulations across the entire 25MHz cellular band have been developed. In an advanced application, a software radio does not just transmit; it characterizes the available transmission channels, probes the propagation path, constructs an appropriate channel modulation, electronically steers its transmit beam in the right direction for systems with antenna arrays, and selects the optimum power level. It does not just receive; it characterizes the energy distribution in the channel and in adjacent channels, recognizes the mode of the incoming transmission, adaptively null interferers, estimates the dynamic properties of multipath propagation, equalizes and decodes the channel codes. The main advantage of software radio is that it supports incremental service enhancements through upgrades to its software. This whole area is not possible without the advancements in DSP technology.
1.5
Brief Overview of Chapters
An overview of the remaining chapters in this book are as follows:
Brief Overview of Chapters

Mobile Unit Speech Video Data
17
ADC DAC
Programmable Processor(s)
Wideband ADC & DAC
RF Conversion
Realtime
Basestation
Service Development Software Offline Programmable Processors Online Programmable Processor(s) Realtime Wideband ADC & DAC RF Conversion
Public Switched Network
Figure 1.8: Software radio architecture Digital signals and systems can either be described as sequences in time or in frequency. In Chapter 2, digital signals are viewed as sequences in time. Digital linear systems are characterized in the time domain by a sequence called the impulse response. We shall discuss the properties of digital signals and systems and their interaction. Signals and systems can also be represented as functions of frequency. Just like the continuous-time Fourier Transform connects the time-domain and frequencydomain view of continuous-time (analog) signals and systems, the discrete-time Fourier Transform (DTFT) is the tool for analyzing discrete-time signals and systems in the frequency domain. Chapter 3 covers the foundamental concepts of the discrete-time Fourier Transform. The DTFT is continuous in frequency. So it is not suitable for computation using a computer. The discrete Fourier transform (DFT) is an approximation of the DTFT that is suitable for numerical computation. The basic characteristics of the DFT and some ways by which the transform can be computed efciently are described in Chapter 4. The application of DFT to linear convolution and spectral analysis is
18 also discussed.
Introduction
Chapter 5 discusses the process of converting a continuous-time signal to a discretetime and discrete-amplitude one and vice versa. Concepts of sampling and quantization and their relation to aliasing are described. The process of reconstructing an analog signal from data sequences is also presented. Special emphasis is placed on the design tradeoffs between sampling rate, quantization resolution, and the complexity of analog lter for sampling and reconstruction. The processing of digital signals are most often performed by digital lters. The different ways by which these FIR and IIR digital lters can be realized by hardware or software will be discussed in Chapter 6. Even though these realizations are based on a single equation, there are surprisingly many variations with different properties. Here the effects of using a nite wordlength to represent numbers and lter coefcients are discussed. The design of the two major types of digital lters: nite impulse response (FIR) and innite impulse response (IIR) lters are coverd in chapters 7and 8. There are many design techniques available. Here we only describe two methods for each type of lters to illustrate the principles involved. Chapters 6 to 8 combined should gives the reader a good understanding of the digital lter approximation and implementation issues. Since this is an introductory course, a number of important but more advanced topics in digital signal processing are not covered. These topics include adaptive ltering, multirate processing, parametric signal modelling and spectral estimation, two (and higher) dimensional digital signal processing, and other efcient fast Fourier transform algorithms. I have also restrict the materials to cover only deterministic signals. However, those interested in exploring these more advanced topics should nd the foundation provided by this book very useful.
Chapter 2 Discrete-time Signals and Systems

2.1 Discrete-time signals
Digital signals are discrete in both time (the independent variable) and amplitude (the dependent variable). Signals that are disrete in time but continuous in amplitude are referred to as discrete-time signals. Discrete-time signals are data sequences. A sequence of data is denoted {x[n]} or simply x [n] when the meaning is clear. The elements of the sequence are called samples. The index n associated with each sample is an integer. If appropriate, the range of n will be specied. Quite often, we are interested in identifying the sample where n = 0. This is done by putting an arrow under that sample. For instance, {x[n]} ={. . . ,0.35,1,1.5,0.6, 2,. . .} The arrow is often omitted if it is clear from the context which sample is x[0]. Sample values can either be real or complex. In the rest of this book, the terms discrete-time signals and sequences are used interchangeably. The time interval between samples is not explicitly shown. It can be assumed to be normalized to 1 unit of time. So the corresponding normalized sampling frequency is 1 Hz. If the actual sampling interval is T seconds, then the sampling frequency is given by 1 fs = T .
20
2 1.8 1.6 1.4 1.2 Amplitude 1 0.8 0.6 0.4 0.2 0 3
Discrete-time Signals and Systems
Index (n)
Figure 2.1: The Unit Impulse Sequence
2.1.1
Some Elementary Sequences
There are some sequence that we shall encounter frequently. They are described here. 2.1.1.1 Unit Impulse Sequence
The unit impulse sequence is dened by [n] = 1, n = 0 0, n = 0 (2.1)
This is depicted graphically in Figure 2.1. Note that while the continuous-time unit impulse function is a mathematical object that cannot be physically realized, the unit impulse sequence can easily be generated. 2.1.1.2 Unit Step Sequence
The unit step sequence is one that has an amplitude of zero for negative indices and an amplitude of one for non-negative indices. It is shown in Figure 2.2.
Discrete-time signals
2 1.8 1.6 1.4 1.2 Amplitude 1 0.8 0.6 0.4 0.2 0 3
21
Index (n)
Figure 2.2: The Unit Step Sequence 2.1.1.3 Sinusoidal Sequences A sinusoidal sequence has the form x[n] = A cos (0 n + ) (2.2)
This function can also be decomposed into its in-phase xi [q] and quadrature xq [n] components. x[n] = A cos cos 0 n A sin sin 0 n = xi [n] + xq [n] This is a common practice in communications signal processing. 2.1.1.4 Complex Exponential Sequences (2.3) (2.4)
Complex exponential sequences are essentially complex sinusoids. x[n] = Aej(o n+) = A cos(0 n + ) + jA sin(0 n + ) (2.5) (2.6)
22
0.5 0.4 0.3 0.2 0.1 Amplitude 0 0.1 0.2 0.3 0.4 0.5
10
15 Index (n)
20
25
30
Figure 2.3: Uniformly distributed random sequence with amplitudes between -0.5 and 0.5. 2.1.1.5 Random Sequences
The sample values of a random sequence are randomly drawn from a certain probability distribution. They are also called stochastic sequences. The two most common distributions are the Gaussian (normal) distribution and the uniform distribution. The zero-mean Gaussian distribution is often used to model noise. Figure 2.3 and gure 2.4 show examples of uniformly distributed and Gaussian distributed random sequences respectively.
2.1.2 Types of Sequences

The discrete-time signals that we encounter can be classied in a number of ways. Some basic classications that are of interest to us are described below.
Real vs. Complex Signals A sequence is considered complex at least one sample is complex-valued.
2.5
23
1.5
Amplitude
0.5
0.5
1.5
10
15 Index (n)
20
25
30
Figure 2.4: Gaussian distributed random sequence with zero mean and unit variance. Finite vs. Innite Length Signals Finite length sequences are dened only for a range of indices, say N1 to N2 . The length of this nite length sequence is given by |N2 N1 + 1|. Causal Signals A sequence x[n] is a causal sequence if x[n] = 0 for n < 0.
Symmetric Signals First consider a real-valued sequence {x[n]}. Even symmetry implies that x[n] = x[n] and for odd symmetry x[n] = x[n] for all n. Any real-valued sequence can be decomposed into odd and even parts so that x[n] = xe [n] + xo [n] where the even part is given by xe [n] = and the odd part is given by xo [n] = 1 (x[n] x[n]) 2 1 (x[n] + x[n]) 2
24
A complex-valued sequence is conjugate symmetric if x[n] = x [n]. The sequence has conjugate anti-symmetry if x[n] = x [n]. Analogous to real-valued sequences, any complex-valued sequence can be decomposed into its conjugate symmetric and conjugate anti-symmetric parts: x[n] = xcs [n] + xca [n] 1 xcs [n] = (x[n] + x [n]) 2 1 xca [n] = (x[n] x [n]) 2 (2.7) (2.8) (2.9)
Periodic Signals A discrete-time sequence is periodic with a period of N samples if x[n] = x[n + kN ] (2.10)
for all integer values of k. Note that N has to be a positive integer. If a sequence is not periodic, it is aperiodic or non-periodic. We know that continuous-time sinusoids are periodic. For instance, the continuous-time signal x(t) = cos (0 t) (2.11) has a frequency of 0 radians per second or f0 = 0 /2 Hz. The period of this sinusoidal signal is T = 1/f0 seconds. Now consider a discrete-time sequence based on a sinusoid. x[n] = cos (0 n) If this sequence is periodic with a period of N samples, then cos 0 (n + N ) = cos 0 n So 0 N must be an integer multiple of 2, or if r is any integer, 0 N =2r 2f0 N =2r r f0 = N (2.13) (2.12)
(2.14)
Since both r and N are integers, a discrete-time sinusoidal sequence is periodic if its frequency is a rational number. Otherwise, it is non-periodic.
Discrete-time signals Example 2.1. Is x[n] = cos

n 8
25 periodic? If so, what is the period?
The sequence can be expressed as x[n] = cos 2 1 16 n
So in this case, f0 = 1/16 is a rational number and the sinusoidal sequence is periodic with a period N = 16. It is interesting to note that for discrete-time sinusoidal sequences, a small change in frequency can lead to a large change in period. For example, a certain sequence has frequency f1 = 31/60. So its period is 60 samples. Another sinusoidal sequence with frequency f2 = 30/60 has a period of only 2 samples since f2 can be simplied to 1/2. Another point to note is that discrete-time sinusoidal sequences with frequencies separated by an integer multiple of 2 are identical. Energy and Power Signals The energy of a nite length sequence x[n] is dened as
N2
E=
n=N1
|x[n]|2
(2.15)
while that for an innite sequence is
E=
n=
|x[n]|2
(2.16)
Note that the energy of an innite length sequence may not be innite. A signal with nite energy is usually referred to as an energy signal. Example 2.2. Find the energy of the innite length sequence x[n] = 2n , n 0 0, n<0
According to the denition, the energy is given by

E=
n=0
2n
=
n=0
1 4
26 To evaluate the nite sum, rst consider

N 1
SN =
n=0
an = 1 + a + a2 + + aN 1
Multiplying this equation by a, we have aSN = a + a2 + + aN and the difference between these two equations give SN aSN = (1 a) SN = 1 aN Hence if a = 1, SN = 1 aN 1a (2.17)
For a = 1, it is obvious that SN = N . For a < 1, the innite sum is therefore S = 1 1a (2.18)
Making use of this equation, the energy of the signal is E= 4 1 = 1 1/4 3
Equations (2.17) and (2.18) are very useful and we shall be making use of them later. The average power of a periodic sequence with a period of N samples is dened as 1 Px = N
N 1
|x[n]|2
n=0
and for non-periodic sequences, it is dened in terms of the following limit if it exists: 1 Px = lim |x[n]|2 K 2K + 1 n=K A signal with nite average power is called a power signal.
K
Discrete-time signals Example 2.3. Find the average power of the unit step sequence u[n]. The unit step sequence is non-periodic, therefore the average power is P = 1 K 2K + 1 lim
27
u2 [n]
n=0
K +1 = lim K 2K + 1 1 = 2 Therefore the unit step sequence is a power signal. Note that its energy is innite and so it is not an energy signal. Bounded Signals A sequence is bounded if every sample of the sequence has a magnitude which is less than or equal to a nite positive value. That is, |x[n]| Bx < Summable Signals A sequence is absolutely summable if the sum of the absolute value of all its samples is nite.
|x[n]| <
n=
A sequence is square summable if the sum of the magnitude squared of all its samples is nite. |x[n]|2 <
n=
2.1.3
Some Basic Operations on Sequences
Scaling Scaling is the multiplication of a scalar constant with each and every sample value of the sequence. This operation is depicted schematically in Figure 2.5. Addition Addition of two sequences usually refers to point-by-point addition as shown in Figure 2.6.
28
x[n]
Ax[n]
Figure 2.5: Scalar Multiplication by A.

x[n] y[n]=x[n]+w[n]
w[n]
Figure 2.6: Point-by-point Addition Delay A unit delay shifts a sequence by one unit of time as shown in Figure 2.7. A sequence x[n] is delayed by N samples to produce y[n] if y[n] = x[n N ].
Up/Down Sampling Down-sampling by a factor of L (a positive integer) is an operation by which only one every L-th sample of the original sequence is kept, with the rest discarded. The schematic for a down-sampler is shown in Figure 2.8. The down-sampled signal y[n] is given by y[n] = x[nM ] Up-sampling is the opposite operation to down-sampling. It increases the number of samples of the original sequence by a certain factor L (a positive integer). This is done by inserting L 1 zeros between each pair of original samples. Figure shows the schematic
x[n]
x[n1]
Figure 2.7: A Unit Delay
29
x[n]
y[n]
Figure 2.8: Down-sampling by a factor of L.

x[n]
y[n]
Figure 2.9: Up-sample by a factor of L. diagram of an up-sampler and the up-sampled sequence y[n] is given by y[n] = x 0,
n L
, n = 0, L, 2L, . . . otherwise
An interpolated signal can be obtained by passing the up-sampled signal through a lowpass lter with an appropriate bandwidth. This process will be discussed in more detail in a later chapter. Modulation Given two sequences x[n] and w[n], and y[n] = x[n] w[n] (2.19)
Then we say that y[n] is w[n] modulated by x[n]. This is analogous to carrier modulation in communication systems. Correlation The correlation, or more precisely cross-correlation, between two nite length data sequences x[n] and w[n] is dened by 1 r= N
N
x[n]w[n]
n=0
(2.20)
if each sequence is of length N . The correlation coefcient r is often used as a measure of how similar the two sequences are. If they are very different, then the value of r is low. The matched lter used in digital communication receivers for optimal detection is also effectively a correlator between the incoming and the template signals.
30
2.2
Discrete-time Systems
A discrete-time system is one that processes a discrete-time input sequence to produce a discrete-time output sequence. There are many different kinds of such systems.
2.2.1 Classication of Systems

Discrete-time systems, like continuous-time systems, can be classied in a variety of ways. Linearity A linear system is one which obeys the superposition principle. For a certain system, let the outputs corresponding to inputs x1 [n] and x2 [n] are y1 [n] and y2 [n] respectively. Now if the input is given by x[n] = Ax1 [n] + Bx2 [n] (2.21)
where A and B are arbitrary constants, then the system is linear if its corresponding output is y[n] = Ay1 [n] + By2 [n] Superposition is a very nice property which makes analysis much simpler. Although many real systems are not entirely linear throughout its operating region (for instance, the bipolar transistor), they can be considered approximately linear for certain input ranges. Linearization is a very useful approach to analyzing nonlinear systems. Almost all the discrete-time systems considered in this book are linear systems. Example 2.4. Are the down-sampler and up-sampler linear systems? Consider the down-sampler y[n] = x[nM ] For input x1 [n], the corresponding output is y1 [n] = x1 [nM ]. For input x2 [n], the output is y2 [n] = x2 [nM ]. Let x[n] be a linear combination of these two input with arbitrary constants A and B so that x[n] = Ax1 [n] + Bx2 [n] The output is given by y[n] = Ax1 [nM ] + Bx2 [nM ] = Ay1 [n] + By2 [n]
Discrete-time Systems Therefore the down-sampler is a linear system. Now consider the up-sampler y[n] = x[n/L], n = 0, L, 2L, . . . 0, otherwise
31
Let y1 [n] and y2 [n] be the outputs for inputs x1 [n] and x2 [n] respectively. For x[n] = Ax1 [n] + Bx2 [n] then the output is y[n] = Ax1 [n/L] + Bx2 [n/L], n = 0, L, 2L, . . . 0, otherwise = Ay1 [n] + By2 [n]
Hence the up-sampler is also linear.
Shift Invariance A shift (or time) invariant system is one that does not change with time. Let a system response to an input x[n] be y[n]. If the input is now shifted by n0 (an integer) samples, x1 [n] = x[n n0 ] (2.22) then the system is shift invariant if its response to x1 [n] is y1 [n] = y[n n0 ] (2.23)
In the remainder of this book, we shall use the terms linear time-invariant (LTI) and linear shift-invariant interchangeably. Example 2.5. A system has input-output relationship given by
n
y[n] =
k=
x[k]
Is this system shift-invariant?
32
If the input is now x1 [n] = x[n n0 ], then the corresponding output is

n
y1 [n] =
k= n
x1 [k] x[k n0 ]
k= nn0
=
k=
x[k]
= y[n n0 ] Therefore the system is shift invariant. Example 2.6. Is the down-sampler a shift invariant system? Let M (a positive integer) be the down-sampling ratio. So for an input x[n] the output is y[n] = x[nM ]. Now if x1 [n] is x[n] delayed by n0 samples, then x1 [n] = x[n n0 ] and the corresponding output is y1 [n] = x[nM n0 ] n0 M = x n M If the system is shift invariant, one would expect the output to be y[n n0 ] = x [(n n0 ) M ] Since this is not the case, the down-sampler must be shift variant. Causality The response of a causal system at any time depends only on the input at the current and past instants, not on any future samples. In other words, the output sample y[n0 ] for any n0 only depends on x[n] for n n0 . Example 2.7. Determine if the following system is causal:
y[n] =
k=
(n k) u[n k]x[k]
33
Note that u[n k] = 0 for n < k because the unit step sequence is zero for negative indices. In other words, for a certain n, u[n k] = 0 for k > n. So the output can be written as
n
y[n] =
k=
(n k) x[k]
So y[n] depends on x[k] for k n and therefore the system is causal.
Stability There are two common criteria for system stability. They are exponential stability and bounded-input bounded-output (BIBO) stability. The rst criterion is more stringent. It requires the response of the system to decay exponentially fast for a nite duration input. The second one merely requires that the output be a bounded sequence if the input is a bounded sequence. Example 2.8. Determine if the system with the following input-output relationship is BIBO stable.
y[n] =
k=
(n k) u[n k]x[k]
Consider input x[n] = [n]. Then
y[n] =
k=
(n k) u[n k][k]
= nu[n] which is unbounded as it grows linearly with n. Therefore the system is not BIBO stable. Example 2.9. Determine if the following system is BIBO stable. Note that this system is an averager, taking the average of the past M samples. 1 y[n] = M
M 1
x[n k]
k=0
34
Let the input |x[n]| B for some nite value B. Consider the magnitude of the output |y[n]| = 1 M 1 M
M 1
x[n k]
k=0 M 1
|x[n k]|
k=0
1 (M B) M = B Hence the output is bounded and the system is BIBO stable. Lossy or Lossless For a passive system that does not generate any energy internally, the output should have at most the same energy as the input. So

|y[n]|
n= n=
|x[n]|2 <
(2.24)
A lossless system is one which the equality holds.
2.2.2
Linear Shift-Invariant Systems
An discrete-time LTI system, like its continuous-time counterpart, is completely characterized by its impulse response. In other words, the impulse response tells us everything we need to know about an LTI system as far as signal processing is concerned. The impulse response is simply the observed system output when the input is an impulse sequence. For continuous-time systems, the impulse function is purely a mathematical entity. However, for discrete-time systems, since we are dealing with sequences of numbers, the impulse sequence can realistically (and easily) be generated. 2.2.2.1 Linear Convolution
Let us consider a discrete-time LTI system with impulse reponse h[n] as shown in Figure 2.10. What would be the output y[n] of the system if the input x[n] is as shown in Figure 2.11?
35
0.8
0.7
0.6
0.5
h[n]
0.4
0.3
0.2
0.1
4 Index (n)
Figure 2.10: Impulse Response of the System Considered in Section 2.2.2.1
2.5
2 Input Signal x[n]
1.5
0.5
4 Index (n)
Figure 2.11: Input Signal Sequence
36
Since the system is linear and time invariant, we can make use of the superposition principle to compute the output. The input sequence is composed of three impulses. Mathematically, it can be expressed as x[n] = [n] + 0.5[n 1] + 2[n 2] Let x1 [n] = [n] x2 [n] = 0.5[n 1] x3 [n] = 2[n 2] and the system response to each of these inputs are respectively y1 [n], y2 [n] and y3 [n]. The sample values of y1 [n] are given by y1 [0] = h[0]x1 [0]=0.8 y1 [1] = h[1]x1 [0]=0.4 y1 [2] = h[2]x1 [0]=0.2 y1 [3] = h[3]x1 [0]=0.1 which is the same as the impulse response since x1 [n] is a unit impulse. Similarly, y2 [1] = h[0]x2 [1]=0.4 y2 [2] = h[1]x2 [1]=0.2 y2 [3] = h[2]x2 [1]=0.1 y2 [4] = h[3]x2 [1]=0.05 and y3 [2] = h[0]x3 [2]=1.6 y3 [3] = h[1]x3 [2]=0.8 y3 [4] = h[2]x3 [2]=0.4 y3 [5] = h[3]x3 [2]=0.2 The system output y[n] in response to input x[n] is therefore, through the superposition principle, given by y[n] = y1 [n] + y2 [n] + y3 [n] = {0.8, 0.8, 2, 1, 0.45, 0.2}
Discrete-time Systems Note that y[0] y[1] y[2] y[3] y[4] y[5] In general, we have y[n] = h[n]x[0] + h[n 1]x[1] + . . . + h[1]x[n 1] + h[0]x[n] or y[n] =
k=0 n
37
= = = = = =
h[0]x[0] h[1]x[0] + h[0]x[1] h[2]x[0] + h[1]x[1] + h[0]x[2] h[3]x[0] + h[2]x[1] + h[1]x[2] h[3]x[1] + h[2]x[2] h[3]x[2]
(2.25)
h[k]x[n k]
n
(2.26)
Alternatively, y[n] =
k=0
x[n]h[n k]
(2.27)
Equations 2.26 and 2.27 are the linear convolution equations for nite length sequences. If the length of h[n] is M and the length of x[n] is N , then the length of y[n] is N +M 1. We can further generalize it for innite length sequences:
y[n] =
k=
h[k]x[n k] x[k]h[n k]
k=
(2.28) (2.29)
These equations are analogous to the linear convolution equation for continuous-time signals. Note that the linear convolution equation comes about as a result of the superposition principles and therefore applies only to LTI systems. The convolution equation for discrete-time signals is also called the convolution sum. It is denoted by . So equations 2.28 and 2.29 can be written as y[n] = x[n] h[n] The convolution sum is one of the most important fundamental equations in DSP.
38 2.2.2.2 Properties of Linear Convolution
The convolution sum has three important properties: 1. Commutative x[n] y[n] = y[n] x[n] 2. Associative (x[n] w[n]) y[n] = x[n] (w[n] y[n]) 3. Distributive x[n] (w[n] + y[n]) = x[n] w[n] + x[n] y[n] 2.2.2.3 Condition for Stability (2.32) (2.31) (2.30)
Since the impulse response completely characterize an LTI system, we should be able to draw conclusions regarding the stability of a system based on its impulse response. We shall consider BIBO stability here. Theorem 2.1. An discrete-time LTI system is BIBO stable if its impulse response is absolutely summable. Proof: Let the input be bounded, i.e. |x[n]| < B < for some nite value B. The magnitude of the output is given by
|y[n]| =
k=
h[k]x[n k] |h[k]| |x[n k|

k=
|h[k]|
k=
So the magnitude of y[n] is bounded if response must be absolutely summable.
k=
|h[k]| is nite. In other words, the impulse
Discrete-time Systems 2.2.2.4 Condition for Causality
39
Theorem 2.2. A discrete-time LTI system is causal if and only if its impulse response is a causal sequence. Proof: Consider an LTI system with impulse response h[k]. Two different inputs x1 [n] and x2 [n] are the same up to a certain point in time, that is x1 [n] = x2 [n] for n n0 for some n0 . The outputs y1 [n] and y2 [n] at n = n0 are given by
y1 [n0 ] =
k= 1
h[k]x1 [n0 k]
=
k=
h[k]x1 [n0 k] +
k=0
h[k]x1 [n0 k]
and
y2 [n0 ] =
k= 1
h[k]x2 [n0 k]
=
k=
h[k]x2 [n0 k] +
k=0
h[k]x2 [n0 k]
Since x1 [n] = x2 [n] for n n0 , if the system is causal, then the outputs y1 [n] and y2 [n] must be the same for n n0 . More specically, y1 [n0 ] = y2 [n0 ]. Now,

h[k]x1 [n0 k] =
k=0 k=0
h[k]x2 [n0 k]
because x1 [n0 k] = x2 [n0 k] for non-negative values of k. Since x1 [n] may not be equal to x2 [n] for n > n0 , we must have
1 1
h[k]x1 [n0 k] =
k= k=
h[k]x2 [n0 k] = 0
which means that h[k] = 0 for k < 0.
40
2.2.3 FIR and IIR Systems

In DSP, nite impulse response (FIR) and innite impulse response (IIR) systems and usually referred to as FIR and IIR lters. An causal FIR lter has impulse response h[n] = {h0 , h1 , h2 , , hM , 0, 0, } where h0 , h1 , h2 , , hM are called lter coefcients. The length of the impulse response is M + 1 and it is called an M -th order lter. Example 2.10. What is the impulse response of a causal LTI system with the following input-output relationship? y[n] = 2x[n] + x[n 1] x[n 3] Impulse response is the systems output when the input is a unit impulse signal. Therefore, h[n] = 2[n] + [n 1] [n 3] = {2, 1, 0, 1} This lter is an FIR lter. The ltering equation is essentially the linear convolution equation. Since the impulse response has a length of M + 1, the output can be computed as (assuming a causal lter):
M
y[n] =
k=0
h[k]x[n k]
(2.33)
The impulse response of an IIR system, however, has innite length. So the ltering equation for a causal lter is given by
y[n] =
k=0
h[k]x[n k]
(2.34)
This equation is not computable in practice due to the innite limit. But there is a type of IIR systems that admits recursive computation of y[n] and so it can be computed. This
41
type of IIR systems are the only IIR lters that are used in DSP. The general input-output relationship of these IIR lters is
N M
y[n] =
k=1
a[k]y[n k] +
k=0
b[k]x[n k]
(2.35)
where {a[k]} and {b[k]} are the lter coefcients. The present output y[n] is a linear combination of previous N outputs and the present and previous M outputs. The order of the lter is M 1 or N , whichever is larger. Example 2.11. A causal LTI system has impulse response h[n] = ah[n 1] + [n] where a is a constant. What is the input-output equation for this system? Since y[n] = h[n] for x[n] = [n], it is obvious from the impulse response that the input-output equation for this system is y[n] = ay[n 1] + x[n] Example 2.12. Find the impulse response of the causal system with the following inputoutput equation: y[n] = 0.25y[n 2] + x[n] The impulse response of this system is given by h[n] = 0.25h[n 2] + [n] and so h[0] = 0.25h[2] + [0]=1 h[1] = 0.25h[1] + [1]=0 h[2] = 0.25h[0] + [2] =(0.25)2 h[3] = 0.25h[1] + [3] =0 . . . . . . Hence h[n] = (0.5)n , n 0, and 0, otherwise even
42
2.3
Summary
In this chapter, we have studied the fundamentals of discrete-time signals and systems in the time domain. Discrete-time signals are essentially a sequence of numbers. Some fundamental sequences, such as the unit impulse, unit step and sinusoidal (both real and complex) sequences are examined because more complex sequences can be expressed as a combination of some or all of these fundamental ones. The most important discrete-time systems are the linear time invariant (LTI) systems. For these systems, the superposition principle applies which leads to the linear convolution sum. This convolution sum is the way by which we derive the output of the system given an input. An LTI system is completely characterized by its impulse response. All the properties of such a system is revealed by its impulse response, including causality and BIBO stability.
Chapter 3 Discrete-time Fourier Analysis

One of the most powerful mathematical tool in signal processing is Fourier analysis. Periodic continuous-time signals can be represented as a Fourier series. For non-periodic signals, the Fourier Transform can be used to provide us with a frequency-domain view of the signals. The Fourier transform has some very nice properties which makes it very useful in theory and in practice. Since in DSP we are dealing with discrete-time signals, there is a need to have a similar set of tools for the frequency analysis of these signals. The Discrete-time Fourier Series (DTFS) and Discrete-time Fourier Transform (DTFT) are such tools. Before we study these techniques in detail, we shall briey review the continuous-time Fourier series and transforms.
3.1
Continuous-time Fourier Analysis
If x(t) is a continuous-time periodic signal, then it can be expressed as a linear combination of weighted exponential functions
x(t) =
c[k]ej2kf0 t
k=
(3.1)
44
Discrete-time Fourier Analysis
where the fundamental period of x(t)is T0 = 1/f0 . The Fourier coefcients c[k] can be computed by T0 1 c[k] = x(t)ej2kf0 t dt (3.2) T0 0 The coefcients c[k] indicate the strength of each frequency component present in the signal. Since x(t) is periodic, the only frequency components present are the harmonics of the fundamental frequency. If x(t) is non-periodic, it can be expressed as
x(t) =
X(f )ej2f t df
(3.3)
where X(f ) =
x(t)ej2f t dt
(3.4)
(3.4) is usually referred to as the forward transform and (3.3) as the inverse transform. X(f ) is the spectrum of the signal x(t). Note that the spectrum of a non-periodic signal is continuous in frequency.
3.2 Discrete-time Fourier Series

Suppose x[n] is a periodic sequence with a period of N samples. It can be expressed as a linear combination of complex exponential functions:
N 1
x[n] =
k=0
c[k]ej2kn/N
n = 0, 1, . . . , N 1
(3.5)
The coefcients c[k] are called discrete-time Fourier coefcients. Note that the complex exponential functions ej2kn/N have a fundamental period of N . In order to obtain an equation for c[k], we start with multiplying both sides of (3.5) by ej2ln/N where l is an integer. So we have
N 1
x[n]e
j2ln/N
=
k=0
c[k]ej2(kl)n/N
(3.6)
Discrete-time Fourier Series Then we sum both sides with respect to n as follows:
N 1 N 1 N 1
45
x[n]e
n=0
j2ln/N
=
n=0 k=0
c[k]ej2(kl)n/N
(3.7)
Since
N 1
an =
n=0
N,
1aN , 1a
a=1 a=1
(3.8)
we can conclude that

N 1
ej2(kl)n/N =
n=0
N, (k l) = 0, N, 2N, . . . 0, otherwise
(3.9)
Therefore the right hand side of (3.7) becomes

N 1 N 1
c[k]
k=0 n=0
ej2(kl)n/N = N c[l]
(3.10)
Put this back into (3.7) and re-arranging, we have 1 N

N 1
c[l] =
x[n]ej2ln/N
n=0
l = 0, 1, . . . , N 1
(3.11)
Note also that these discrete-time Fourier coefcients are also periodic with a period of N. Example 3.1. Determine the discrete-time Fourier coefcients of the following periodic sequence with a period of 4 samples: x[n] = {. . . , 0, 1, 1, 0, . . .} The coefcients are calculated using (3.2) with N = 4. 1 c[k] = 4
3
x[n]ej2kn/4
n=0
1 = x[0] + x[1]ejk/2 + x[2]ejk + x[3]ej3k/2 4 1 = x[1]ejk/2 + x[2]ejk 4
46 So c[0] c[1] c[2] c[3] = = = =
0.5 0.25 (1 j) 0 0.25 (1 + j)
3.2.1 Average Power

The average power of a periodic sequence x[n] is dened as 1 P = N where N is the period of x[n]. We want to express the average power in terms of c[k] instead of x[n]. 1 P = N = 1 N
N 1 N 1
|x[n]|2
n=0
(3.12)
|x[n]|2
n=0 N 1
x[n]x [n]
n=0 N 1 N 1
1 = N
x[n]
n=0 k=0
c [k]ej2kn/N
Interchanging the order of summation, we have

N 1
P =
k=0 N 1
c [k]
1 N
N 1
x[n]ej2kn/N
n=0
=
k=0 N 1
c [k]c[k] |c[k]|2
k=0
(3.13)
Discrete-time Fourier Series
47
Now we have established that the average power can be calculated using either the signal sequence or its discrete-time Fourier coefcients. The relationship 1 N is called the Parsevals relation. Example 3.2. Verify the Parsevals relation for the periodic sequence in Example 3.1. Solution: Using the sequence sample values, we have
3 N 1 N 1
|x[n]| =
n=0 k=0
|c[k]|2
(3.14)
P =
1 4
|x[n]|2
n=0
1 = (1 + 1) 4 1 = 2 Now using the Fourier coefcients, we have

3
P =
k=0
|c[k]|2
= (0.5)2 + |0.25 (1 j)|2 + |0.25 (1 + j)|2 = 0.25 + 0.125 + 0.125 = 0.5 Thus the Parsevals relation is veried.
3.2.2
Normalized Frequency
Each Fourier coefcient is associated with a particular frequency. From (3.5), we can see that c[k] corresponds to the complex sinusoid ej2kn/N = cos (2kn/N ) + jsin (2kn/N ) = cos (k n) + jsin (k n) (3.15) (3.16)
48 with an angular frequency k = For k = 0, 1, . . . , N 1, 0 k < 2 2k N
(3.17)
So the Fourier series expansion only gives us frequency components within a 2 range. k is the normalized (angular) frequency because we implicity assume that the time spacing T between samples of the data sequence is 1 second. If the actual value of T is known, then we can obtain the de-normalized frequencies k = and the range given by 0 k < 2 T (3.19) 2k NT (3.18)
3.2.3
Power Density Spectrum
Since each Fourier coefcient c[k] is associated with a particular frequency, |c[k]|2 is proportional to the power of that particular frequency component within the discretetime signal. A plot of |c[k]|2 against k shows the distribution of power as a function of frequency. The sequence |c[k]|2 , k = 0, 1, . . . , N 1 is called the power density spectrum of the signal. The power density spectrum provide us with a frequency-domain view, as opposed to the time-domain view, of the signal. Example 3.3. Plot the power density spectrum of the periodic signal given in Example 3.1. The power density spectrum is given by
|c[k]|2 = The plot is shown in Figure 3.1.
0.25, 0.125, 0, 0.125
Discrete-time Fourier Transform

Power Density Spectrum
49
0.25 Magnitude Squared Fourier Coefficients
0.2
0.15
0.1
0.05
0 1
0.5
0.5
1.5 k
2.5
3.5
Figure 3.1: Power Density Spectrum in Example 3.3.
3.3
Discrete-time Fourier Transform
Using the discrete-time Fourier series, we can analyze periodic discrete-time signals in the frequency domain. In order to deal with the more general non-periodic signals, we need to use the discrete-time Fourier Transform (DTFT). It is dened as
X () =
n=
x[n]ejn
(3.20)
Notice that the DTFT is continuous in the variable , while x[n] is discrete-time. Also note that in order to compute the DTFT, we need to know the whole data sequence (n = to ).
50
3.3.1 Convergence
The DTFT is a sum of an innite sequence. As we know, not all innite sums converge. Therefore we would like to establish the conditions under which the DTFT exists. The DTFT is said to exist if the nite sum
K
XK () =
n=K
x[n]ejn
(3.21)
converges to X () in some sense as K approaches innity. 3.3.1.1 Uniform Convergence
XK () as dened in (3.21) is said to converge uniformly if

K
lim |X () XK ()| = 0
(3.22)
Notice that in this case we require that the absolute difference between X () and XK () becomes smaller and smaller as K increases at every single frequency . Consider
|X ()| =
n=
x[n]ejn |x[n]|
n=
(3.23) (3.24)
So the magnitude of X () is nite if x[n] is absolutely summable. 3.3.1.2 Mean Square Convergence
If DTFT is only applicable to absolutely summable sequences, then its use will be very limited. Most of the data sequences in practice are not absolutely summable, but they are square summable. By relaxing the convergence criteria, we can extend the DTFT to include square summable sequences.
Discrete-time Fourier Transform The mean square convergence criterion is dened as

K
51
lim
|X () XK ()|2 d = 0
(3.25)
Notice that this criterion requires that the sum of the differences squared become smaller as K increases. In other words, the energy of the error decreases when more terms are included. However, it does not mean that for any particular frequency, the difference is reduced. This fact is illustrated through the following example. Example 3.4. Consider the following discrete time signal x[n] = The DTFT of x[n] is given by
sin c n , n c ,
n=0 n=0
X () =
n=
x[n]ejn 1, || c 0, c < ||
= x[n] is square summable since
|x[n]|
n=
1 = 2 1 = 2 c =
|X ()|2 d
c
d
c
However, it is not absolutely summable. So its DTFT exists only in the mean square sense. In Figure 3.2 we plotted XK () for several different values of K to illustrate mean square convergence. Here we use c = 0.5. The DTFT of x[n] is an ideal low-pass lter. As K increases, it can be observed that the overall error becomes smaller due to smaller errors away from the jump discontinuity at c . However, the peak of the ripple in the passband remains very much the same height regardless of K. This effect is due to mean square convergence at a jump discontinuity.
52
1.4
1.2
K=21 K=11
K=51
Magnitude
0.8
0.6
0.4
0.2
0.5
1.5 2 Normalized frequency
2.5
3.5
Figure 3.2: Finite term approximations to the DTFT. At the end of the 19-th century, Albert Abraham Michelson conducted an experiment to measure the speed of light using an interferometer. This experiment became known as the Michelson-Morley experiment and Michelson subsequently became the rst American to win the Nobel prize. He also built a device that performs Fourier analysis on light. He observed that where there is a jump discontinuity in the signal, some sort of ringing occurs. He could not explain the phenomenon and consulted Josiah Willard Gibbs. Gibbs eventually came up with an explanation and published it in 1899 in Nature. So this became known as the Gibbs phenomenon.
3.3.2
Energy Density Spectrum
The energy density spectrum of a non-periodic discrete-time signal is given by |X ()|2 . It tells us the energy distribution of the signal in frequency.
Discrete-time Fourier Transform The magnitude and phase spectra are given by |X ()| and X () respectively.
53
The energy of a nite-energy sequence can be calculated using the data sequence itself or from the energy density spectrum.
E =
n=
|x[n]|2
(3.26) (3.27)
1 = 2 This is the Parsevals relation for DTFT.
|X ()|2 d
3.3.3
3.3.3.1
Properties
Periodicity
The DTFT is periodic with a period of 2. This can easily be veried as follows. Let k be any integer,
X ( + 2k) =
n=
x[n]ej(+2k)n x[n]ejn ej2kn

n=
= =
n=
x[n]ej()n (3.28)
= X ()
This means that the DTFT of a discrete-time sequence is only unique within a (normalized) frequency range of 2. The DTFT is repeated once every 2. Figure 3.3 shows the periodic repetition of the magnitude spectrum of discrete-time sequences. 3.3.3.2 Symmetry
If x[n] is a real-valued discrete-time signal, then its DTFT has conjugate symmetry, i.e. X () = X () (3.29)
54
Replica of |X(w)|
|X(w)|
Replica of |X(w)|
....
3 2 2 3
....
Figure 3.3: Periodic nature of the Fourier spectrum of discrete-time sequences. or X () = X () (3.30) Furthermore, if x[n] is real-valued and even, then its DTFT is also real-valued and even. If x[n] is real-valued and odd, then its spectrum is imaginary-valued and odd. 3.3.3.3 Linearity DTFT is a linear transform since it obeys the superposition principle. For any scalar constants a and b, and data sequences x[n] and y[n] where x[n] X () y[n] Y () are DTFT pairs, ax[n] + by[n] aX () + bY () 3.3.3.4 Time Shifting (3.33) (3.31) (3.32)
If the data sequence x[n] is shifted by an integer index of n0 , the effect in the frequency domain is a shift in phase by an amount n0 at frequency . x[n n0 ] ejn0 X () (3.34)
Discrete-time Fourier Transform 3.3.3.5 Time Reversal
55
Reversing the time sequence has the effect of reversing the DTFT. x[n] X () 3.3.3.6 Frequency Shifting (3.35)
Analogous to time-shifting, if the DTFT X () is shifted in frequency by a constant 0 , then the effect in the time domain is a multiplication by a complex sinusoid. ej0 n x[n] X ( 0 ) (3.36)
A particular case of frequency shifting is when the time sequence is modulated by a sinsoidal signal. x[n]cos (0 n) 1 {X ( + 0 ) + X ( 0 )} 2 (3.37)
This is commonly encountered in communication systems. The effect of modulation is a shift of the original DTFT up and down the frequency axis by an amount equal to the modulating frequency. 3.3.3.7 Linear Convolution
One of the most useful property of DTFT is converting a linear convolution in the time domain into multiplication in the frequency domain. Let x1 [n] X1 () x2 [n] X2 () Consider the linear convolution of x1 [n] and x2 [n]: y[n] = x1 [n] x2 [n]
(3.38) (3.39)
=
k=
x1 [k]x2 [n k]
(3.40)
56
Multiplying both sides by ejn and sum over all n, we have

y[n]e
n=
jn
=
n= k=
x1 [k]x2 [n k] ejn
(3.41)
The left hand side of the above equation is the DTFT of y[n] which we shall denote by Y (). Therefore,

Y () =
k=
x1 [k]
n=
x2 [n k]ejn x2 [m]ej(m+k)
m= jk m=
=
k=
x1 [k] x1 [k]e
k=
x2 [m]ejm (3.42)
= X1 () X2 ()
3.3.4 Inverse DTFT

The time domain sequence can be recovered given its DTFT. Comparing (3.43) with (3.1), we realize that (3.43) can be viewed as a (continuous) Fourier series representation of X (). So the inverse DTFT should take the form of (3.2). The inverse DTFT is given by x[n] = 1 2
X () ejn d
(3.43)
It can be shown that this is indeed the inverse transform by substitution of the above expression into (3.20).
3.3.5
Autocorrelation Sequence
In many applications, the autocorrelation sequence of a discrete-time signal x[n] has to be computed. It is dened as
rxx [l] =
n=
x[n]x[n l]
l = 0, 1, 2, . . .
(3.44)
Summary where the parameter l is called the lag. Note that rxx [0] is the energy of the signal. We can express the autocorrelation with a lag of l as a convolution sum.
57
rxx [l] =
n=
x [n] x [ (l n)]
= x[n] x[n] If x[n] is a real-valued sequence and the DTFT of x[n] is X (), then the DTFT of x[n] is X (). So DT F T {x[n] x[n]} = X () X () = |X ()|2 Hence we can compute the energy density spectrum of real-valued signals by taking the DTFT of the autocorrelation sequence of the signal. This is particularly useful for random signals.
3.4
Summary
Signals and systems can be viewed from the frequency domain perspective. We showed that a periodic time sequence can be expressed as a sum of a nite number of sinuoids (the Fourier series) similar to that of continuous-time signals. Non-periodic discretetime signals can also be expressed as the sum of an innite number of sinusoids (the discrete-time Fourier Transform). The main difference between continuous-time Fourier Transforms (CTFT) and DTFT is that the DTFT is periodic. Properties of the DTFT are otherwise similar to those for CTFT. Two of the more important properties are linear convolution and the Parsevals theorem.
58
Chapter 4 Discrete Fourier Transform and Its Computation

A discrete-time sequence can be converted to its frequency domain representation by the Discrete-time Fourier Transform. The frequency variable is a continuous variable. While it is not difcult to compute the DTFT at a particular frequency, it is impossible to compute the whole DTFT numerically using a computer. One way to approximate the DTFT is to sample it at regular frequency intervals similar to the way we sample time domain signals. We call the sampled version of the DTFT the Discrete Fourier Transform (DFT). In this chapter, we are concerned about the DFT and its implementation and applications to computing the linear convolution and to spectral analysis.
4.1
Discrete Fourier Transform
4.1.1 Sampling the DTFT

Suppose we take N evenly spaced frequency samples from the DTFT X () of a discretetime sequence x [n]. The sampled values are given by X [k] = X ()|=2k/N
=
n=
x [n] ej2kn/N
(4.1)
60
m=1 N
Discrete Fourier Transform and Its Computation

m=2 m=3
Samples with same exponential factor are added up
Figure 4.1: Illustration of Equation 4.3 Using 3 Data Blocks. for 0 k N 1. We can divide the summation in (4.1) into blocks of N terms:
(m+1)N 1
X [k] =
m= n=mN N 1
x [n] ej2kn/N x [n mN ] ej2kn/N

n=0 m=
x[N1+2N]
x[N1+N]
x[0+2N]
x[1+2N]
x[N1]
x[0+N]
x[1+N]
x[0]
x[1]
(4.2)
(4.3)
(4.3) is obtained from (4.2) by interchanging the inner and outer summations with appropriate change of variables. The term in bracket in (4.3) is the sum of all x [n] terms with the same exponential factor. This equation is illustrated in Figure with 3 blocks. Suppose x [n] is of nite length and the length of x [n] is N , then
x [n mN ] = x [n]
m=
(4.4)
since x [n] = 0 for n < 0 and n N . So (4.3) becomes

N 1
X [k] =
n=0
x [n] ej2kn/N
(4.5)
61
for k = 0, 1, . . . , N 1. This is called the N -point Discrete Fourier Transform of x [n]. X [k] are known as the DFT coefcients. Note that the DFT transforms a length-N data sequence into a length-N sequence of DFT coefcients. In other words, if we need more sample points in the frequency domain, we need more data points in the time domain.
4.1.2 Inverse DFT

Note that in (4.3), the term m= x [n mN ] (denoted as xp [n]) is periodic with a period of N . Therefore it can be expressed as a discreter Fourier series, i.e.
N 1
xp [n] =
k=0
c [k] ej2kn/N
where c [k] are the Fourier series coefcients. They are given by 1 c [k] = N
N 1
xp [n] ej2kn/N
n=0
1 = X [k] N using (3.11) in Section 3.2 and (4.5). If x [n] has a nite length of N samples, then x [n] = xp [n] and we have 1 x [n] = N
N 1
0nN 1
X [k] ej2kn/N
k=0
(4.6)
for 0 n N . (4.6) allows us to recover the time sequence x [n] from the DFT coefcients X [k]. It is therefore called the inverse Discrete Fourier Transform (IDFT). The IDFT can be obtained through a forward DFT. Taking the complex conjugate of (4.6) we have 1 x [n] = N =
N 1
X [k] ej2kn/N
k=0
(4.7) (4.8)
1 N
DF T (X [k])
62
So by computing the DFT of the complex conjugate of the DFT coefcients, we can obtain the complex conjugate of x [n] and hence x [n] itself.
4.1.3 Matrix Form

The factor ej2/N is central to the DFT equations. It is usually denoted as WN = ej2/N and so the DFT equation (4.5) becomes
N 1
(4.9)
X [k] =
n=0
kn x [n] WN
(4.10)
for 0 k N 1. The inverse DFT equation (4.6) then becomes 1 x [n] = N for 0 n N 1. (4.10) can be expressed in matrix form: X = WN x Or X [0] X [1] X [2] . . . X [N 1] 1 1 1 . . . 1 1 WN 2 WN . . . 1 2 WN 4 WN . . .
2(N 1) N 1 kn X [k] WN k=0
(4.11)
(4.12)
.. .
1
N WN 1 2(N 1) WN . . . (N 1)(N 1)
x [0] x [1] x [2] . . . x [N 1]
N 1 WN 1 WN
WN
(4.13) Notice that some of the entries in the WN matrix are the same. For instance, it can easily be shown that k+N/2 k (4.14) WN = WN Similarly for (4.11),
1 x = WN X
(4.15)
Discrete Fourier Transform where it can be shown that

1 WN =
63
1 W N N
(4.16)
Example 4.1. Derive the W4 matrix and hence compute the 4-point DFT of x [n] = 1, 2, 3, 4
1 Then derive the W4 matrix and recover x [n] from its DFT coefcients.
Solution:
0 W4 W0 W4 = 4 0 W4 0 W4 0 W4 W0 W4 = 4 0 W4 0 W4 k+N k Since WN = WN and WN
0 W4 1 W4 2 W4 3 W4 0 W4 1 W4 2 W4 3 W4
0 W4 2 W4 0 W4 6 W4 0 W4 2 W4 0 W4 6 W4
0 W4 1 1 1 1 3 1 j 1 j W4 = 2 W4 1 1 1 1 9 1 j 1 j W4 0 W4 1 1 1 1 3 1 j 1 j W4 = 2 W4 1 1 1 1 9 1 j 1 j W4
k+N/2
k = WN , 6 W4 9 W4 2 W4 3 W4
= = = =
2 W4 1 W4 0 W4 1 W4
and 1 1 1 1 1 1 j 1 j 2 X = 1 1 1 1 3 1 j 1 j 4 10 2 + 2j = 2 2 2j
64
1 Now WN = 1 WN N
Discrete Fourier Transform and Its Computation and so

1 W4
1 1 1 1 1 1 j 1 j = 4 1 1 1 1 1 j 1 j
To recover x [n], 1 1 1 1 10 1 j 1 j 2 + 2j 1 x = 2 4 1 1 1 1 1 j 1 j 2 2j 1 2 = 3 4
4.1.4
Zero-padding
If we have a length-L data sequence x [n] and want to perform an N -point DFT where L < N , we can append x [n] with D = N L zeros to make it of length N . x [n] = {x0 , x1 , , xL1 } xD [n] = {x0 , x1 , , xL1 , 0, . . . , 0} (4.17) (4.18)
Note that by zero-padding a signal, we do not alter its DTFT, i.e. XD () = X ():
N 1
XD () =
n=0 L1
xD [n] ejn
L+D1
(4.19) +
n=L
=
n=0 L1
xD [n] e
jn
xD [n] ejn
(4.20)
=
n=0
x [n] ejn
(4.21) (4.22)
= X () since x [n] = xD [n] for n = 0, 1, . . . , L 1 and xD [n] = 0 for n L. However, the DFT coefcients XD [k] and X [k] are not the same.
65
4.1.5 Circular Shift

Consider a length N data sequence x [n] which is zero outside the range 0 n N 1 . This sequence is then time shifted by n0 samples so it becomes x [n n0 ] for some non-zero values of n0 . What would be the N -point DFT of x [n n0 ]? Lets use an example. x [n] = {1, 2, 3, 4} and n0 = 2. So the time-shift sequence is xs [n] = {0, 0, 1, 2, 3, 4} (4.24) (4.23)
and we want a 4-point DFT of xs [n]. From (4.3), the DFT of xs [n] is the same as the DFT of length-4 sequence xc [n] = {3, 4, 1, 2} (4.25) Note that xc [n] is a circularly shifted version of x [n]: xc [n] = x [ n n0 where n
N N]
(4.26)
denotes n modulo N . In other words, xc [n] = x [n n0 ] , n0 n N 1 x [N n0 + n] , 0 n < n0 (4.27)
So a time-shift is equivalent to a circular shift as far as the DFT is concerned.
4.1.6 Properties
In this section, some of the important properties of the DFT are summaried. Periodicity Consider X [k + N ] which is given by
N 1
X [k + N ] =
n=0 N 1
x [n] ej2kn/N ej2N n/N x [n] ej2kn/N

n=0
(4.28) (4.29) (4.30)
= X [k]
66
This implies that X [k] is periodic with a period of N , even though x [n] is aperiodic. In general, X [k] = X [k + pN ] p = 1, 2, . . . (4.31) Thus for a discrete-time signal, we cannot obtain the full spectrum (frequencies from negative innity to innity). Linearity The DFT is a linear transformation. Hence the superpostion principle holds. If both x1 [n] and x2 [n] are of the same length and x1 [n] X1 [k] x2 [n] X2 [k] then x [n] = ax1 [n] + bx2 [n] X [k] = aX1 [k] + bX2 [k] where a and b are arbitrary constants. Parsevals Relation The Parsevals relation holds for DFT, similar to the continuoustime Fourier Transform and the DTFT. So the energy of the signal can be calculated from the time sequence or from the DFT coefcients.
N 1
|x(n)|2 =
n=0
1 N
N 1
|X(k)|2
k=0
Time Shifting We noted in Section 4.1.5 that time delaying a sequence is equivalent to circularly shift a sequence as far as the DFT is concerned. The DFT of the circularly shifted sequence is given by x[ n m
N] km WN X [k]
(4.32)
where X [k] are the DFT coefcients of the original non-shifted sequence. In other words, a time (circular) shift in the time domain leads to a linear phase shift of an amount ej2km/N in the frequency domain. Frequency Shifting A circular shift in the DFT coefcients will have similar effect time domain: mn WN x [n] X [ k m N ] (4.33)
67
Symmetry If x [n] is real-valued, then its DFT has the following symmetry properties: {X [k]} = {X [k]} {X [k]} = {X [k]} Furthermore, if x [n] is real and even, then X [k] is also real and even; and if x [n] is real and odd, then X [k] is imaginary and odd. If x [n] is purely imaginary, then {X [k]} = {X [k]} {X [k]} = {X [k]} Also, if x [n] is imaginary and even, then X [k] is also imaginary and even; and if x [n] is imaginary and odd, then X [k] is real and odd. The proof is straight-forward and is left as an exercise to the reader. (4.36) (4.37) (4.34) (4.35)
4.1.7
Circular Convolution
Let X1 [k] and X2 [k] be N -point DFTs of sequences x1 [n] and x2 [n] respectively so that
N 1
X1 [k] =
n=0 N 1
x1 [n] ej2kn/N x2 [n] ej2kn/N

n=0
(4.38) (4.39)
X2 [k] =
68
Now consider Y [k] = X1 [k] X2 [k] for k = 0, 1, . . . , N 1. Inverse DFT gives 1 y [m] = N 1 = N = Now 1 N
N 1
X1 [k] X2 [k] ej2km/N

k=0 N 1 N 1 N 1
(4.40) x2 [l] ej2kl/N

l=0
x1 [n] e
k=0 N 1 n=0 N 1
j2kn/N
ej2km/N (4.41)
N 1
x1 [n]
n=0 l=0
x2 [l]
k=0
ej2k(mnl)/N
(4.42)
N 1
ej2k(mnl)/N =
k=0
N,
1aN , 1a
a=1 a=1
(4.43)
where a = ej2(mnl)/N N, l = m n + pN = m n = 0, otherwise So y [n] =

n=0 N 1
(4.44)
N
(4.45)
x1 [n] x2 [ m n
N]
(4.46)
The circular convolution y [n] between two length-N sequences x1 [n] and x2 [n] is dened by (4.46). An N -point circular convolution is denoted by y [n] = x1 [n] x2 [n] (4.47)
Therefore a circular convolution can be computed by multiplying the DFTs of the two data sequences followed by an inverse DFT of the product. Example 4.2. Compute the circular convolution of the following length-4 sequences using DFT. x [n] = {1, 2, 1, 2} h [n] = {1, 2, 3, 4} Solution:
Discrete Fourier Transform Circular convolution using DFT: X [k] = DF T {x [n]} = {6, 0, 2, 0} H [k] = DF T {h [n]} = {10, 2 + j2, 2, 2 j2} and Y [k] = H [k] X [k] = {60, 0, 4, 0} So the circular convolution of x [n] and h [n] gives y [n] = IDF T {Y [k]} = {16, 14, 16, 14}
69
The circular convolution operation in (4.46) can also be expressed in matrix form. x1 [0] y [0] x2 [0] x2 [N 1] x2 [1] y [1] x2 [1] x2 [0] x2 [2] x1 [1] = (4.48) . . . . . .. . . . . . . . . . . . y [N 1] x2 [N 1] x2 [N 2] x2 [0] x1 [N 1] Example 4.3. Check the answer obtained in Example 4.2 by using the denition (4.46). Solution: Using the matrix form, 1 2 y [n] = 3 4
4 1 2 3
3 4 1 2
2 1 3 2 4 1 1 2
1+8+3+4 16 2 + 2 + 4 + 6 14 = 3 + 4 + 1 + 8 = 16 4+6+2+2 14
which is the same as the answer obtained in Example 4.2.
4.1.8 Linear Convolution Using DFT

Circular convolution does not have any real meaning, unlike linear convolution. Therefore we are interested in the relationship between linear and circular convolution. If we can obtain the linear convolution through circular convolution, then we can make use of DFT to perform linear convolution calculations. Suppose two sequences x [n] and h [n] have lengths N and M respectively. The linear convolution of these two sequences will yield a result which is a sequence y [n] of length
70
(N + M 1). So if we zero-pad x [n] and h [n] to length (N + M 1) and perform a (N + M 1)-point circular convolution, the result should be the same as that obtained by linear convolution. Example 4.4. Compute the linear convolution of the data sequences x [n] = {1, 2, 0, 1} h [n] = {2, 2, 1, 1} Solution: The linear convolution result should be of length (4 + 4 1) = 7. Therefore 3 zeros are need for each sequence. xD [n] = {1, 2, 0, 1, 0, 0, 0} hD [n] = {2, 2, 1, 1, 0, 0, 0} So the circular convolution of xD [n] and hD [n] is given by 2 0 0 0 1 1 2 1 2 2 2 0 0 0 1 1 2 2 + 4 1 2 2 0 0 0 1 0 1 + 4 y [n] = 1 1 2 2 0 0 0 1 = 1 + 2 + 2 0 1 1 2 2 0 0 0 2 + 2 0 0 1 1 2 2 0 0 1 0 0 0 1 1 2 2 0 1 It can be veried that this is the result of linear convolution. Overlap-Add Method In many cases we have a relatively short sequence h [n] (e.g. an impulse response sequence) and we need to compute the output y [n] (e.g. the output of a linear system) given a very long sequence x [n] (e.g. the input to the linear system). If we used the method above, the output will only be obtained when the whole input is collected and the convolution sum computed. If x [n] is indeed very long, the delay will be long as well. If x [n] is innitely long, then the method cannot be applied. An alternative way is to compute the linear convolution through circular convolution is to divide the input into nite length blocks. The idea is to obtain a block of output samples at a time.
2 6 5 5 4 1 1
71
Given a long sequence x [n], divide it into non-overlapping sub-sequences xr [n] of length L so that x [n] =
r=0
xr [n rL]
(4.49)
Let h [n] be a shorter sequence of length M . Append (L 1) zeros to h [n] and (M 1) zeros to each sub-sequence xr [n] so that they are both of length N = M + L 1. Compute the linear convolution of xr [n] and h [n] (using DFT if desired), yielding an output yr [n] of length N . y [n] = x [n] h [n]
(4.50) h [n] (4.51) (4.52)
=
r=0
xr [n rL] yr [n rL]
r=0
= where
yr [n] = xr [n] h [n]
(4.53)
Notice in (4.52) that each block of yr [n] is of length N while the shift between yr [n] and yr1 [n] is L samples. Therefore there is an overlap of (N L) or (M 1) samples between subsequent blocks of yr [n]. According to (4.52), these overlapped sample values should be added to give the correct y [n]. This process is shown in Figure 4.2 for three blocks of data. Example 4.5. Compute the linear convolution of the following sequences using 4-point DFTs. h [n] = {3, 2, 1} x [n] = {1, 2, 3, 4, 4, 3, 2, 1} Solution: Length of the shorter sequence is M = 3. Use 4-point FFTs, therefore N = 4. So L=N M +1=43+1=2 Hence x [n] is to be divided into blocks of 2 samples.
72
L x[n] M h[n] L y 0[n] M1 x 0[n] L x [n]
1

L x [n]
2
L x m [n]
y [n]
1
y [n]
2
y[n]
Figure 4.2: Illustration of the Overlap-add Method. We have h [n] = {3, 2, 1, 0} x0 [n] = {1, 2, 0, 0} x1 [n] = {3, 4, 0, 0} x2 [n] = {4, 3, 0, 0} x3 [n] = {2, 1, 0, 0} and so Y0 [k] = H [k] X0 [k] Y1 [k] = H [k] X1 [k] Y2 [k] = H [k] X2 [k] Y3 [k] = H [k] X3 [k] = = = = {18, 2 6j, 2, 2 + 6j} {42, 2 14j, 2, 2 + 14j} {42, 2 14j, 2, 2 + 14j} {18, 2 6j, 2, 2 + 6j} H [k] = {6, 2 2j, 2, 2 + 2j} X0 [k] = {3, 1 2j, 1, 1 + 2j} X1 [k] = {7, 3 4j, 1, 3 + 4j} X2 [k] = {7, 4 3j, 1, 4 + 3j} X3 [k] = {3, 2 j, 1, 2 + j}
M1
M1
Discrete Fourier Transform Inverse FFTs give us y0 [n] y1 [n] y2 [n] y3 [n] Align them properly and add: 2 18 11 4 12 17 10 3 + 6 7 4 1 3 8 14 20 23 21 16 10 4 1 which gives us y [n] = {3, 8, 14, 20, 23, 21, 16, 10, 4, 1}. 3 8 5 9 = = = = {3, 8, 5, 2} {9, 18, 11, 4} {12, 17, 10, 3} {6, 7, 4, 1}
73
Overlap-Save Method There is an alternative way to the overlap-add method. Consider a length-3 sequence h [n] and a length-4 sequence x [n]. the result of the linear convolution of h [n] and x [n], which is of length 6, is given by yL [0] yL [1] yL [2] yL [3] yL [4] yL [5] = = = = = = h [0] x [0] h [0] x [1] + h [1] x [0] h [0] x [2] + h [1] x [1] + h [2] x [0] h [0] x [3] + h [1] x [2] + h [2] x [1] h [1] x [3] + h [2] x [2] h [2] x [3] (4.54) (4.55) (4.56) (4.57) (4.58) (4.59)
Compare this with the circular convolution of these sequences (after appending a zero to h [n]): yc [0] yc [1] yc [2] yc [3] = = = = h [0] x [0] + h [1] x [3] + h [2] x [2] h [0] x [1] + h [1] x [0] + h [2] x [3] h [0] x [2] + h [1] x [1] + h [2] x [0] h [0] x [3] + h [1] x [2] + h [2] x [1] (4.60) (4.61) (4.62) (4.63)
74 Obviously,
yc [0] = yL [0] yc [1] = yL [1] However, yc [2] = yL [2] yc [3] = yL [3]
(4.64) (4.65) (4.66) (4.67)
In general, if the length of the shorter sequence is M and we use an N point circular convolution, then the last L = (N M + 1) samples are the same as the linear convolution output. In other words, the rst (M 1) samples of the circular convolution output are useless as far as the linear convolution is concerned. Therefore they can be discarded and only the remaining output samples are kept. All these saved output samples can now be concatenated to form the correct linear convolution output. This process is illustrated in Figure 4.3 for three data blocks. This method is called Overlap-save. Example 4.6. Repeat example 4.5 using the overlap-save method. Solution: Again L = 2. But instead of appending zeros to each block, we have x0 [n] = {0, 0, 1, 2} x1 [n] = {1, 2, 3, 4} x2 [n] = {3, 4, 4, 3} x3 [n] = {4, 3, 2, 1} x4 [n] = {2, 1, 0, 0} X0 [k] = {3, 1 + 2j, 1, 1 2j} X1 [k] = {10, 2 + 2j, 2, 2 2j} X2 [k] = {14, 1 j, 0, 1 + j} X3 [k] = {10, 2 2j, 2, 2 + 2j} X4 [k] = {3, 2 j, 1, 2 + j} y0 [n] = {5, 2, 3, 8} y1 [n] = {14, 12, 14, 20} y2 [n] = {19, 21, 23, 21} y3 [n] = {16, 18, 16, 10} y4 [n] = {6, 7, 4, 1}
Each Xi [k] , i = 0, . . . , 4 is multiplied with H [k] to give Y0 [k] = {18, 2 + 6j, 2, 2 6j} Y1 [k] = {60, 8j, 4, 8j} Y2 [k] = {84, 4, 0, 4} Y3 [k] = {60, 8j, 4, 8j} Y4 [k] = {18, 2 6j, 2, 2 + 6j}
Retaining only the last L = 2 samples of each yi [n] gives us the linear convolution output y [n] = {3, 8, 14, 20, 23, 21, 16, 10, 4, 1} which is the same as that obtained using overlap-add.

M1 L x 0[n] L x [n]
1
75
L x [n]
2 m
L x [n]
x[n]
h[n] M1 L y 0[n] M1 y [n]

1
Discard
y[n]
Figure 4.3: Illustration of Overlap-save Method.
4.1.9
Computational Complexity
Consider the DFT equation:

N 1
X [k] =
n=0
To compute each X [k], we need to perform N complex multiplications (the multiplication of two complex numbers) and N 1 complex additions (the addition of two complex numbers). So the complete N -point DFT requires N 2 complex multiplications and N (N 1) complex additions in general. Note that the multiplication of two complex numbers (a + jb) (c + jd) = (ac bd) + j (bc + ad)
M1 y [n]
2
x [n] ej2kn/N
(4.68)
76
involves 4 real number multiplications and 2 real number additions. The addition of two complex numbers involves 2 real number additions. Therefore an N -point DFT involves 4N 2 real multiplications and 2N (N 1) + 2N 2 real additions. To compute an N point circular convolution using the equation
N 1
y [n] =
n=0
x1 [n] x2 [ m n
N]
(4.69)
we require N 2 multiplications and N (N 1) additions. If the sequences x1 [n] and x2 [n] are real-valued, then these multiplications and additions involve only real numbers. If the circular convolution is computed using N -point DFTs, then it involves two DFTs to transform x1 [n] and x2 [n] to X1 [k] and X2 [k], and one inverse DFT to transform Y [k] = X1 [k] X2 [k] to the result y [n]. Since the IDFT equation is essentially the same as the DFT equation, their computation complexities are the same. So the two DFTs and one IDFT together requires 3N 2 complex multiplications and 3N (N 1) complex additions. To compute Y [k], we need another N complex multiplications. So the total number of complex multiplications is 3N 2 + N and the number of complex additions is 3N (N 1). Note that in general, the DFT coefcients are complex-valued even though the data sequences x1 [n] and x2 [n] are real-valued.
4.2
The Fast Fourier Transform
Based on the computational complexity analysis in Section 4.1.9, the DFT does not seem very useful for performing circular convolution (and hence linear convolution using overlap-add and overlap-save). In the early 1960s, Cooley and Tukey developed an efcient algorithm to compute N -point DFTs where N = 2r for positive integers r. It is called the fast Fourier transform (FFT). This algorithm has a signicant impact on the use of DFT in various computations and essentially launched the eld of digital signal processing.
77
4.2.1 Radix-2 Decimation-in-time Algorithm

The basic idea of Cooley and Tukeys FFT algorithm, which is now known as the Radix-2 decimation-in-time algorithm, is to rewrite the DFT equation (4.5) into two parts:
N/21 N/21 2nk x [2n] WN n=0
X [k] =
k WN n=0
2nk x [2n + 1] WN
(4.70)
for k = 0, 1, . . . , N 1 and N = 2r for some positive integer r. The rst part of (4.70) is the DFT of the even sequence and the second part, is the DFT of the odd sequence. Notice that the factor W 2nk appears in both DFTs and need only be computed once. Also note that
2nk nk WN = WN/2
(4.71)
Let x1 [n] = x [2n] and x2 [n] = x [2n + 1], and x1 [n] X1 [k] x2 [n] X2 [k] (4.70) can be written as
k X [k] = X1 [k] + WN X2 [k]
(4.72)
Or
k X [k] = X1 [k] + WN X2 [k] N k X k+ = X1 [k] WN X2 [k] 2
(4.73) (4.74)
for k = 0, 1, . . . , N 1 since WN 2 twiddle factor.
k+N/2
k k = WN . The complex factor WN is known as the
Each of the two subsequences x1 [n] and x2 [n] can be further broken down into even and odd sequences. This process can continue until only 2-point DFTs are left. So each N/2point DFT is obtained by combining two N/4-point DFTs, each of which is obtained by combining two N/8-point DFTs, etc. There are a total of r stages since N = 2r . Figure 4.4 shows the concept for an 8-point FFT. The algorithm is called decimation-in-time because we break up the original data (or time) sequence into two reduced number of samples (decimation). Radix-2 refers to the fact that 2-point DFTs are the basic computational blocks of the algorithm.
78
x[0] x[4]

2point DFT
X[0] X[1] X[2]
Combine 2point
x[2] x[6]
2point DFT
DFTs Combine 4point
X[3] X[4]
x[1] x[5]
2point DFT
DFTs Combine 2point

X[5] X[6] X[7]
x[3] x[7]
2point DFT
DFTs
Stage 1
Stage 2
Stage 3
Figure 4.4: The three stages of an 8-point Radix-2 decimation-in-time FFT.
Buttery Operation
Computation of two-point DFTs are trivial. It involves
0 X2p [0] = x2p [0] + W2 x2p [1] = x2p [0] + x2p [1] 1 X2p [1] = x2p [0] + W2 x2p [1] = x2p [0] x2p [1]
(4.75) (4.76)
Combining the N -point DFT coefcients to form an N -point DFT involves the compu2 tation of (4.73) and (4.74). The basic operation is illustrated in Figure 4.5. It is called a buttery operation and is a basic building block of the FFT structure.
Example 4.7. Derive the Radix-2 decimation-in-time FFT structure for a 4-point DFT. Solution:
79
a+W Nb
b WN
k
aWN b
Figure 4.5: Decimation-in-time FFT Buttery Operation The 4-point DFT is rst decomposed:
3
X [k] =
n=0
nk x [n] W4
k = 0, 1, 2, 3
= =
0 2k k 3k x [0] W4 + x [2] W4 + x [1] W4 + x [3] W4 0 2k k 0 2k x [0] W4 + x [2] W4 +W4 x [1] W4 + x [3] W4 2pt DF T 2pt DF T
k = X1 [k] + W4 X2 [k]
k = 0, 1, 2, 3
More explicitly, the second stage of the FFT is governed by the equations:
k X [k] = X1 [k] + W4 X2 [k] , k X [k + 2] = X1 [k] W4 X2 [k] ,
k = 0, 1 k = 0, 1
The rst stage consists of two 2-point DFTs governed by the equations:
0 X1 [0] = x [0] + W4 x[2] 2 0 X1 [1] = x [0] + W4 x [2] = x [0] W4 x [2]
and
0 X2 [0] = x [1] + W4 x[3] 2 0 X2 [1] = x [1] + W4 x [3] = x [1] W4 x [3]
80
x[0]

X[0]
x[2]
W4 x[1]
X[1]
W4 x[3] W4
0
X[2]
W4
X[3]
Figure 4.6: Signal Flow Graph of the Radix-2 Decimation-in-time 4-point FFT Algorithm. Figure 4.7: Signal Flow Graph of the Radix-2 Decimation-in-time 8-point FFT Algorithm. Signal Flow Graph The FFT structure is usually represented in a signal ow graph. The signal ow graph for the radix-2 DIT FFT algorithm in Example 4.7 is shown in Figure 4.6. In the graph, the arrows represent the direction of the signal ow, the nodes represent adders and the multiplication factors are written on the segment where it is to be performed. Figure 4.7 shows the signal ow graph of an 8-point radix-2 DIT FFT. Bit Reversal Note that if we want the DFT coefcients to come out in the natural order, the input sequence has to be rearranged as shown in Figures 4.6 and 4.7. This reordering is known as bit reversal. We can see why it is called bit reversal if we represent the index sequence in binary form. This is illustrated by Table. A simple algorithm for generating a list of bit reversed numbers is attributed to Bruneman. It goes as follows: 1. Start with {0, 1}. Multiply by 2 to get {0, 2}.
The Fast Fourier Transform 2. Add 1 to the list of numbers obtained above. In this case, it is {1, 3}. 3. Append the list in step 2 to that in step 1 to get {0, 2, 1, 3}.
81
4. The list obtained in step 3 now becomes the starting list in step 1. The steps are repeated until the desired length of the list is obtained. It should be pointed out that more efcient algorithms have been developed for bit reversal. But we shall not cover them here. Computational Complexity An N -point FFT consists of N/2 butteries per stage with log2 N stages. Each buttery has one complex multiplication and two complex additions. Thus there are a total of (N/2)log2 N complex multiplications compared with N 2 for DFT, and N log2 N complex additions compared with N (N 1) for the DFT. This represents a substantial saving when N is large.
4.2.2
Radix-2 Decimation-in-frequency Algorithm
Another radix-2 FFT algorithm can be derived by partitioning the data sequence into two halves as follows:
N/21 N 1 kn x [n] WN n=0
X [k] =
+
n=N/2
kn x [n] WN
(4.77)
instead of odd and even sequences. The second term on the right of (4.77) can be expressed as
N 1 kn x [n] WN n=N/2 N/21
=
n=0 N/21
x n+
N k(n+N/2) WN 2 N kN/2 kn WN WN 2 N kn WN 2
(4.78)
=
n=0
x n+
N/21 k n=0
(4.79)
= (1)
x n+
(4.80)
82 So (4.77) becomes
N/21
N/21 kn x [n] WN
X [k] =
n=0 N/21
+ (1)
k n=0
x n+ N 2
N kn WN 2 (4.81)
=
n=0
x [n] + (1)k x n +
kn WN
The FFT coefcient sequence can be broken up into even and odd sequences and they have the form
N/21
X [2k] =
n=0 N/21
x [n] + x n +
N 2 N 2
2kn WN
=
n=0 N/21
x [n] + x n +
kn g1 [n] WN/2 n=0 N/21
kn WN/2
(4.82) N 2 N 2
(2k+1)n
X [2k + 1] =
n=0 N/21
x [n] x n +
WN
=
n=0 N/21
x [n] x n +
n WN n=0 N/21
2kn n WN WN
x [n] x n +
N 2
kn WN/2
=
n=0
kn g2 [n] WN/2
(4.83)
for k = 0, 1, . . . , N/2 1 where g1 [n] = x [n] + x n + g2 [n] = N 2 N x [n] x n + 2 (4.84)

n WN
(4.85)
83
The data sequences involved in the even and odd coefcient sequences X [k] and X [2k + 1] can again be divided into the rst and second halves and the same procedure repeated until only 2-point DFTs are required. This algorithm is called decimation-in-frequency (DIF) FFT because it is the DFT coefcients that are divided into odd and even halves. Figure shows the three stages of an 8-point radix-2 decimation-in-frequency FFT. Note that the data sequence appears in its natural order whereas the FFT output occurs in bit reversed order. The computational complexity is the same as the decimation-in-time algorithm.
Buttery Operation The computation of the g1(n) and g2(n) sequences involve the buttery operation as shown in Figure 4-5. It is similar to the buttery for the decimationin-time FFT except for the position of the twiddle factor.
a a+b
WN b 1 W N (ab)
k
Figure 4.8: Radix-2 Decimation-in-Frequency Buttery Example 4.8. Derive the 4-point radix-2 DIF FFT algorithm structure.
Signal Flow Graph The even and odd coefcient sequences can each be divided into the rst and second halves and the same procedure repeated until only 2-point DFTs are required. Figure 4-6 shows the three stages of an 8-point radix-2 decimation-infrequency FFT. Note that the data sequence appears in its natural order whereas the FFT output occurs in bit reversed order. The computational complexity is the same as the decimation-in-time algorithm.
84
4.2.3 Other Fast Algorithms

There are literally tens of different variations of the basic radix-2 FFT algorithms. Some further computational saving can be obtained by using radix-4 or mixed-radix (splitradix) variants. Other algorithms deal with cases where the number of data points is not a power of two, typically when N is a prime number. One of the most well known in this category is the Winograd Fourier transform (WFT). It does not make use of the buttery and the idea is completely different from the Cooley and Tukey algorithm. It is very fast in the sense that it requires fewer multiplications and is particularly useful for short length transforms. But both the mathematics involved and the implementation are considerably more complex. Thus they are outside of the scope for this introductory course.
4.3
Spectral Analysis
The two main applications of DFT are in the computation of the output of a linear system through the linear convolution and the analysis of the spectral content of a signal. In this section we shall concentrate on the latter application. Recall that the DFT coefcients are the sampled values of the DTFT of a given discretetime signal. X [k] = X ()|=2k/N (4.86) for k = 0, 1, . . . , N 1. For a sampling frequency fs , each index k corresponds to the actual frequency k k = fs (4.87) N
4.3.1
Spectral Leakage
Let us consider two examples. Example 4.9. Compute the DFT of a length-32 sinusoidal sequence of frequency 10 Hz, sampled at 64 Hz and plot the magnitude spectrum. Solution: Since it is a single sinusoid, we expect the DFT coefcients to have non-zero values only at N k= fs
Spectral Analysis
85
1 0.8 0.6 0.4 0.2 Amplitude 0 0.2 0.4 0.6 0.8 1
10 Hz sinusoid, sampled at 64 Hz
0.05
0.1
0.15
0.2
0.25 time (s)
0.3
0.35
0.4
0.45
0.5
(a) 10 Hz sinusoid sampled at 64 Hz
20 18 16 14 12 10 8 6 4 2 0
Magnitude of DFT coefficients
10
15 DFT index k
20
25
30
(b) Magnitude spectrum of the sampled 10 Hz sinusoid with N=64
Figure 4.9: Signal and Spectrum in Example 4.9.
86
according to (4.87). Now, N = 32, fs = 64 and = 10. So this gives k = 5. The sampled signal and its magnitude spectrum is shown in Figure 4.9. There is another non-zero value at k = 32 5 = 27 due to the symmetry of the DFT of a real-valued signal. Example 4.10. Compute the DFT of a length-32 sinusoidal sequence of frequency 11 Hz, sampled at 64 Hz and plot the magnitude spectrum. Solution: Again, use k= N fs
we get k = 11/2 = 5.5. The peak is expected to be in between two DFT coefcients. The magnitude spectrum is plotted in Figure 4.10. We may expect the magnitude spectrum in example 4.10 to consist of only two major components similar to that in example 4.9. However, this is not the case. The energy contained in a single frequency is now spread to many frequencies, with the total energy remaining the same (by Parsevals theorem). This effect is called energy leakage. Examining the data sequences reveals why leakage occurs. The signal in example 4.9, shown in Figure 4.9, begins and terminates at the beginning and end of a cycle respectively. However, the signal in example 4.10, shown in Figure 4.10, ends mid-way through a cycle. Therefore the signal that is being transformed has an abrupt jump to zero. This is because the assumptions we made in obtaining the DFT equation is that the signal is zero outside the N samples (see Section 4.1.1). This jump in amplitude contributes to the extra frequency components in the spectrum as shown. In order to avoid leakage, we need to choose N , the size of the DFT, to be contain an integer multiple of cycles of the sinusoid. While we can do that for a single sinusoid, when we have a signal that contains many sinusoidal components, leakage becomes unavoidable.
4.3.2
Rectangular Windowing
Another way of looking at a nite length sinusoidal signal g [n] is that it is a produce of an innite length sinusoid x [n] of the same frequency multiplied by a nite length
Spectral Analysis
87
1 0.8 0.6 0.4 0.2 Amplitude 0 0.2 0.4 0.6 0.8 1
11 Hz sinusoid, sampled at 64 Hz
0.05
0.1
0.15
0.2
0.25 time (s)
0.3
0.35
0.4
0.45
0.5
(a) 11 Hz sinusoid sampled at 64 Hz
12
10
Magnitude of DFT coefficients
10
15 DFT index k
20
25
30
(b) Magnitude Spectrum of sampled sinusoid with N=64
Figure 4.10: Signal and Spectrum in Example 4.10.
88 rectangular window w [n]. Here
w [n] =
1, N1 n N2 0, otherwise
(4.88)
For convenience, N1 is usually assumed to be zero and N2 = N 1 for a window of length N . The DTFT of w [n]is
N 1
W () =
n=0
w [n] ejn
1 ejN 1 ej ejN/2 ejN/2 ejN/2 = ej/2 (ej/2 ej/2 ) sin (N/2) j(N 1)/2 = e sin (/2) = So the magnitude spectrum is given by |W ()| = Since sin for very small , |W ()| = N/2 =N /2 sin (N/2) sin (/2)
(4.89)
(4.90)
(4.91)
for close to zero. This is shown in Figure 4.11. It consists of one main lobe and many sidelobes. The spectrum is zero when the numerator in (4.90) is zero, i.e. sin (N/2) = 0 This occurs when = for k = 1, 2, . . . . 2k N (4.92)
Spectral Analysis
89
The main lobe is of height N in agreement with (4.91). The base width of the main lobe is 4/N . The width of the main lobe may also be dened by the 3 dB width where |W ()| drops by half. In this case, = or 2fs N (4.93)
fs N when it is denormalized by the sampling frequency fs . f =

12 Magnitude Spectrum of Rectangular Window of Length 10
(4.94)
10 Main lobe
8 Magnitude of DTFT
Sidelobes
0 1
0.8
0.6
0.4
0.2 0 0.2 Normalized frequency (Hz)
0.4
0.6
0.8
Figure 4.11: Magnitude Spectrum of a Rectangular Window of length 10.
4.3.3
Frequency Resolution
Suppose we wish to analyze the spectrum of a signal that is composed of the sum of two sinusoids. What is the minimum data length N required so that we can see the peaks of these sinusoids distinctly?
90
Two distinct peaks in a DFT spectrum will merge with one another if the main lobes corresponding to the two sinusoids overlap substantially. This happens when the difference in frequency of the sinusoidal components is less than the main lobe width. Using the 3-dB main lobe width denition, the width is given by (4.94). If the desired frequency resolution is fr which may be the frequency separation of the sinusoidal components, then the minimum number of samples has to ensure that fr f where f is the main lobe width. So the number of data samples required is given by fs N (4.95) fr In other words, the higher the resolution, the smaller the value of fr , and so the more data samples we need. Example 4.11. A signal consists of four sinusoidal frequencies: 1, 1.5, 2.5 and 2.75 kHz. It is sampled at 10 kHz. How many samples are needed for the spectrum to exhibit four peaks? Solution: If a rectangular window is used, the main lobe width is given by f = 2.75 2.5 = 0.25kHz since this is the smallest frequency separation between the four frequency components. Let fr = f , then 10 fs N = = 40 f 0.25 So we need at least 40 samples of data. Thus the main lobe width of the DTFT of the windowing function determines the achieveable frequency resolution. The sidelobes determine the amount of frequency leakage. Obviously we want to increase the main lobe width and the amplitudes of the sidelobes. For the rectangular window, the amplitude of the rst (and thus highest) sidelobe compared to the main lobe is given by W (3/N ) 2 W (0) 3 or in terms of dB, R = 20 log10 W (3/N ) 13.46dB W (0) (4.97) (4.96)
Spectral Analysis
91
To obtain better performance, one way is to design a window to give us the desired parameters. Hamming Window The Hamming window of length L is dened as 0.54 0.46 cos 0,
2n L1
w [n] =
, 0nL1 otherwise
(4.98)
A 51 point Hamming window is shown in Figure 4.12.

51point Hamming Window
0.8
Amplitude
0.6
0.4
0.2
10
15
20
25 Index (k)
30
35
40
45
50
Figure 4.12: A Hamming Window of Length 51. The effective width of the main lobe of the Hamming window is given by fw = 2fs L (4.99)
This translates to a main lobe width of 8/L. So for the same data length, the Hamming window can provide only half the frequency resolution compared to a rectangular window. The advantage of the Hamming window is in the reduction of spectral leakage. The ratio of the rst side lobe peak amplitude to the main lobe is R = 40dB.
92
There exists other windows such as the Hanning window, Blackman window and the Kaiser window. They mostly possess the raised-cosine shape as the Hamming window. In general, they give better values of R but worse frequency resolution compared to the rectangular window.
4.3.4
Computational Frequency Resolution
An estimation of the spectrum of a discrete-time signal x [n] can be obtained through a segment of data. The data segment is usually obtained by windowing x [n] to obtain g [n]. The estimate is the DTFT of g [n]. In practice, DFT is used which gives us a sampled version of the DTFT. Since FFT is much more computationally efcient, it is generally used to compute the DFT. However, the length of the data segment being analyzed may not be of the convenient length for FFT. This problem can be overcome by appending g [n] with the necessary number of zeros. Suppose the length of g [n] is L and it is appended with zeros to produce a sequence gD [n] of length N where L < N . Performing an L-point DFT on g [n] with give us L DFT coefcients. Performing an N -point FFT on gD [n] will give us N FFT coefcients. Since N > L,does that mean we can increase the frequency resolution by appending zeros? The answer is a denite no. As discussed in Section 4.3.3, the frequency resolution is determined by the data length of g [n]. This resolution will not change by appending zeros to g [n]. The apparent increase in the number of frequency points provided by the FFT of gD [n] is actually a re-sampling of the L frequency samples. No extra information can be obtained through re-sampling. This situation is similar to that of a digital camera. The maximum resolution of the image that can be acquired is determined by the resolution of the image sensor on the camera. Any increase in resolution by interpolation is purely articial. No extra details of the image could be obtained through pure interpolation alone. Example 4.12. A signal is composed of two sinusoids of frequencies 1 Hz and 2 Hz. It is sampled at a rate of 5 Hz. How many samples are required to resolve two frequencies using DFT? Solution: If a rectangular window is used, L 5 fs = =5 f (2 1)
Spectral Analysis samples. However, if a Hamming window is used, L samples. 2 (5) 2fs = = 10 fw (2 1)
93
4.3.5
Non-stationary Signals
So far we have only considered signals that are stationary, i.e. their parameters (such as frequency components) do not change over time. However, many signals encountered in real life do change over time. Examples include speech and music signals. They are generally known as non-stationary signals. If we want to know the spectral content of this kind of signals, we cannot simply analyze one segment of the signal. Consider a chirp which has the form x [n] = A cos 0 n2 (4.100)
This signal is shown in Figure 4.13. It is obvious that the instantaneous frequency (the frequency of the sinusoid at a particular time) changes with time. So analyzing a long segment of the chirp signal will show the averaged frequency over that segment which is a misleading result. One way to tackle the spectral analysis of non-stationary signals is to assume that within a reasonably short time, the signal remains stationary. So we divide the signal into sufciently small segments by windowing. The spectra obtained from the spectral analysis of each segment are then displayed side by side to show the changes. To show the slowly varying spectrum with better time resolution, these segments may overlap one another. The resulting display is called a spectrogram. Figure 4.14 shows a spectrogram of the chirp signal using 256-point FFTs, Hamming window, and 250-point overlap. In this gure, the darker the colour, the higher the magnitude. It clearly shows that the instantaneous frequency of the signal increases linearly with time. The length of the windowed segment determines the frequency resolution while the amount of overlap between subsequent segments determines the time resolution. To compute a spectrogram with good frequency and time resolution will require a lot of computation (FFTs). It is common to either relax the frequency resolution or the time resolution requirements in order to lessen the computational requirements. Spectrograms with good frequency resolution are called narrowband spectrograms while those with good time resolution are called wideband spectrograms.
94

1kHz Sampling rate from DC for 0.5s
1.5
0.5
Amplitude
0.5
1.5
0.05
0.1
0.15
0.2
0.25 Time (s)
0.3
0.35
0.4
0.45
0.5
Figure 4.13: A Linear Chirp Signal.
4.4
Summary
In this chapter we discussed ways to compute the spectrum of a discrete-time signal. The Discrete Fourier Transform is introduced as a sampled version of the discrete-time Fourier transform of a discrete-time signal. It can be more efciently computed through fast algorithms called fast Fourier Transforms. The two main applications of DFT/FFT is in performing linear convolution (to obtain the output of discrete-time linear systems) and in the spectral analysis of signals. For the rst application, the overlap-add and overlapsave methods could be employed in cases where the input signal could be innitely long. For spectral analysis, we discussed the relationship between the frequency resolution and the data length required. The analysis of non-stationary signals is also briey discussed.
Summary
95
Spectrogram of linear chirp 500 450 400 350 300 Frequency 250 200 150 100 50 0
0.2
0.4
0.6
0.8 Time
1.2
1.4
1.6
Figure 4.14: Spectrograms of the Chirp Signal.
96
Chapter 5 Sampling, Reconstruction and Quantization

Most of the signals encountered in engineering applications are analog. In order to process analog signals using digital techniques they must rst be converted into digital signals. Digital processing of analog signals proceeds in three stages: 1. The analog signal is digitized. Digitization involves two processes: sampling (digitization in time) and quantization (digitization in amplitude). This whole process is called analog-to-digital (A/D) conversion. 2. The digitized signal is processed by the appropriate DSP algorithms. 3. The results or outputs of the processing are converted back into analog signals through interpolation. This process is called digital-to-analog (D/A) conversion. Figure 5.1 illustrates the main functional blocks involved in the process.
5.1 The Sampling Process

The analog signal is sampled once every T seconds, resulting in a sampled data sequence. The sampler is generally assumed to be ideal in that the value of the signal at an instant (an innitely small time) is taken. A real sampler, of course, cannot achieve that and
98
AnalogtoDigital Conversion
Sampling, Reconstruction and Quantization

Analog Reconstruction
Analog Input
lowpass filter
Sample and Hold
Quantizer
DSP
Digital to Analog Converter
lowpass filter
Analog Output
Antialiasing Prefilter
Antiimage Postfilter
Figure 5.1: Digital Processing of Analog Signals.
Ideal Sampler
Analog Signal x(t)
Sampled Discretetime Signal x(nT)
Figure 5.2: Schematic Diagram of Ideal Sampling
the switch in the sampler (see Figure 5.2) is actually closed for a nite, though very small, amount of time. This is analogous to a camera with a nite shutter speed. Even if a camera can be built with an innitely fast shutter, the amount of light that can reach the lm plane will be very small indeed. However, for simplicity in analysis the sampling process is usually considered to be close enough to the ideal.
It should be pointed out that throughout our discussions we shall assume that the sampling period T is constant. In other words, the spacing between the samples is regular. This is called uniform sampling. Although irregularly sampled signals can, under suitable conditions, be converted to uniformly sampled ones, the concept and mathematics are beyond the scope of this introductory course.
The Sampling Process
99
5.1.1 Spectra of Sampled Signals

The most important parameter in the sampling process is the sampling period T , or the sampling frequency or sampling rate fs which is dened as fs = 1 T (5.1)
Sampling frequency is in units of samples per second or Hertz. The sampled sequence x [n] is given by x [n] = x (nT ) n = . . . , 1, 0, 1, . . . (5.2)
where x (t) is the analog signal. Alternatively, we can view the sampled signal as an analog signal which can be expressed as:
xa (t) = x (t)
n=
(t nT )
(5.3)
where (t) is the impulse function. The impulse train
s (t) =
n=
(t nT )
(5.4)
is a period continuous-time function. So it can be expressed as a Fourier series
s (t) =
m=
ck ej2mt/T
(5.5)
The Fourier coefcents are ck 1 = T 1 = T 1 = T

T
s (t) ej2kt/T dt
0 T
(5.6) (5.7) (5.8)
(t) ej2kt/T dt
0
100
So combining the above equations, we have xa (t) = x (t) 1 T
ej2mt/T
m=
(5.9)
From the frequency shifting property of the Fourier Transform, we know that if we have the Fourier Transform pair x (t) X () (5.10) then we have the pair x (t) ej0 t X ( 0 ) (5.11) for some constant 0 . Therefore, since xa (t) consists of x (t) multiplied by an innite number of ej0 t where 0 = 2m/T , the Fourier transform of xa (t) is given by 1 Xa () = T where
X ( ms )
m=
(5.12)
2 = 2fs (5.13) T is the sampling frequency in radians per second. In other words, the spectrum of the sampled signal consists of frequency shifted replica of the spectrum of the original continuoustime signal. s = Figure 5.3 shows the magnitude spectrum of the sampled signal when the original continuoustime signal is low-pass bandlimited. A low-pass bandlimited signal with a bandwidth of W has a spectrum which is non-zero for the frequency range 0 || W .
5.1.2
Aliasing and The Uniform Sampling Theorem
Each spectrum replica in Figure 5.3 is separated from the adjacent replicas by the sampling frequency. If we reduce the sampling frequency, the replicas will be closer to each other. If the sampling frequency is so low that the replicas overlap each other, as shown in Figure , then the spectrum will be corrupted. We will not be able to reconstruct the original analog signal from the spectrum if it is corrupted. This situation is known as aliasing. Note that this happens when s < 2W .

Magnitude Spectrum X(+2 s ) X(+s ) X( ) X(s ) X(2 s )
101
...
2 s s W 0 W s 2 s
...
Figure 5.3: Spectrum of a Sampled Bandlimited Analog Signal. Uniform Sampling Theorem If an analog signal x (t) is bandlimited to W Hz, then it can be uniquely specied by samples taken uniformly at a rate of at least 2W samples per second, i.e. the sampling frequency and sampling period are given by fs 2W 1 T 2W This is known as the uniform sampling theorem. A sampling rate of 2W for an analog signal of bandwidth W is called the Nyquist rate. A signal is under-sampled if it is sampled below the Nyquist rate; it is critically sampled if it is sampled at the Nyquist rate; it is over-sampled if it is sampled at higher than the Nyquist rate. The sampling theorem provides us with the lower bound on the sampling rate. If the sampling is too sparse, then important information will be missing in the sampled signal. However, if the sampling rate is too high, then a large amount of data will have to be processed within a given time frame. Thus the processing speed of the digital signal processor gives us the upper bound on the sampling rate. It is interesting to note that even though this theorem is usually referred to as Shannons sampling theorem, it was originated by the British mathematicians E.T. and J.M. Whittaker. In the Russian literature, this theorem was introduced to communications theory by Kotelnikov and took its name from him. C.E. Shannon used it to study what is now (5.14) (5.15)
102

Magnitude Spectrum X(+2 s ) X(+s ) X( )
X(s ) X(2 s )
...
Aliasing
...
Figure 5.4: Aliasing of the Frequency Spectrum. known as Information Theory in the 1940s. Therefore in the mathematics and engineeirng literatures sometimes it is also called WKS sampling theorem after Whittaker, Kotelnikov and Shannon. Effect of Aliasing The effect of aliasing on an input signal in the time domain can be demonstrated by sampling a sine wave of frequency fb using different sampling frequencies fs . Figure 5.5 shows such a sinusoidal function sampled at three different rates: fs = 4fb , fs = 2fb , and fs = 1.5fb . In the rst two cases, if we join the sample points using straight lines, it is obvious that the basic up-down nature of the sinusoid is still preserved by the resulting triangular wave as shown in Figure 5.6. If we pass this triangular wave through a low-pass lter, a smooth interpolated function will result. If the low-pass lter has the appropriate cut-off frequency, the original sine wave can be recovered. This will be discussed in detail in Section 5.2. For the last case, the sampling frequency is below the Nyquist rate. We would expect aliasing to occur. This is indeed the case. If we join the sampled points together, it can be observed that the rate at which the resulting function repeats itself differs from the frequency of the original signal. In fact, if we interpolate between the sample points, a smooth sinusoidal function with a lower frequency results, as shown in Figure 5.7. Therefore it is no longer possible to recover the original sine wave from these sampled points. We say that the higher frequency sine wave now has an alias in the lower frequency sine wave inferred from the samples. In other words, these samples are no longer
103
Sampling at 4 times the frequency of the sinusoid
1 0.8 0.6 0.4 0.2 amplitude 0 0.2 0.4 0.6 0.8 1
0.1
0.2
0.3
0.4
0.5 time
0.6
0.7
0.8
0.9
(a) Oversampled 2 times
Sampling at 2 times the frequency of the sinusoid
1 0.8 0.6 0.4 0.2 amplitude 0 0.2 0.4 0.6 0.8 1
0.1
0.2
0.3
0.4
0.5 time
0.6
0.7
0.8
0.9
(b) Critically Sampled
Sampling at 1.5 times the frequency of the sinusoid
1 0.8 0.6 0.4 0.2 amplitude 0 0.2 0.4 0.6 0.8 1
0.1
0.2
0.3
0.4
0.5 time
0.6
0.7
0.8
0.9
(c) Under-sampled
Figure 5.5: A Sinusoid Sampled at 3 Different Rates.
104
1 0.8 0.6 0.4 0.2 amplitude 0 0.2 0.4 0.6 0.8 1
interpolated
0.1
0.2
0.3
0.4
0.5 time
0.6
0.7
0.8
0.9
(a) Over-sampled case
1 0.8 0.6 0.4 0.2 amplitude 0 0.2 0.4 0.6 0.8 1 interpolated
0.1
0.2
0.3
0.4
0.5 time
0.6
0.7
0.8
0.9
(b) Critically Sampled Case
Figure 5.6: Interpolation of Sample Points without Aliasing.

Sampling at 1.5 times the frequency of the sinusoid
105
1 0.8 0.6 0.4 0.2 amplitude 0 0.2 0.4 0.6 0.8 1
1 Hz
2 Hz
0.1
0.2
0.3
0.4
0.5 time
0.6
0.7
0.8
0.9
Figure 5.7: Effect of Aliasing. uniquely representative of the input signal and therefore any subsequent processing will be invalid. The baseband aliased frequency is a frequency within the range f2s , f2s . It is the frequency range by which the samples will be reconstructed to form an analog signal and is known as the Nyquist interval. If a signal is under-sampled, the frequency of the signal lies outside of this frequency range. Since the spectrum of the sampled signal is replicated up and down the frequency axis, a replica appears in the baseband aliased frequency range. This frequency will be the one that is reconstructed. Using the above example, the sinusoid has a frequency of 2 Hz and it is sampled at 3 Hz. So the Nyquist interval is [1.5, 1.5]. The 2 Hz is reected into this range as 1 Hz shown in Figure. In general, the baseband aliased frequency fa of any undersampled signal of frequency f at a sampling rate fs is given by fa = f nfs (5.16)
106
Sampling, Reconstruction and Quantization f2s ,

fs 2
where n is an integer so that fa lies within
Example 5.1. A 7 Hz sinewave is sampled at a rate of 6 Hz. What will be the baseband aliased frequency of this sinewave? Solution:
f = 7 fs = 6 So the apparent frequency is f = f fs = 1 Hz Example 5.2. Determine the Nyquist rate for the signal: x (t) = 3 cos (t) + cos (3t) Suppose this signal is sampled at half the Nyquist rate, determine the signal xa (t) that will be aliased with x (t). Solution: The signal has two frequency components: f1 = 0.5 Hz f2 = 1.5 Hz So the Nyquist sampling rate, which is twice the highest frequency, is 3 Hz. Now if fs = 1.5 Hz, the Nyquist interval is [0.75, 0.75]. f1 falls within the Nyquist interval and is not aliased. But f2 will be aliased with the frequency f2a = f2 fs = 1.5 1.5 = 0 Hz So the signal that will be aliased with x (t) is xa (t) = 3 cos (t) + cos (2f2a t) = 1 + 3 cos (t)
107
5.1.3 Anti-aliasing Prelters
The uniform sampling theorem assumes that the signal is strictly bandlimited. In the real world, typical signals have a wide spectrum and are not bandlimited in the strict sense. For instance, we may assume that 20kHz is the highest frequency the human ears can detect. Thus we want to sample at a frequency slightly above 40kHz (say, 44.1kHz as in compact discs) as dictated by the Sampling Theorem. However, the actual audio signals normally have a much wider bandwidth than 20kHz. We can ensure that the signal is bandlimited at 20kHz by low-pass ltering. This low-pass lter is usually called antialias lter. Anti-aliasing lters are analog lters as they process the signal before it is sampled. In most cases they are also low-pass lters unless bandpass sampling techniques are used. (Bandpass sampling is beyond the scope of this course.) An ideal low-pass lter has a at passband and the cut-off is very sharp. Since the cut-off frequency of this lter is half of that of the sampling frequency, the resulting replicated spectrum of the sampled signal do not overlap each other. Thus no aliasing occurs. Practical low-pass lters cannot achieve the ideal characteristics. What are the implications? Firstly, this would mean that we have to sample the ltered signals at a rate that is higher than the Nyquist rate to compensate for the transition band of the lter. The bandwidth of a low-pass lter is usually dened as the 3-dB point (the frequency at which the magnitude response is 3dB below the peak level in the passband). But signal levels below 3dB are still quite signicant for most applications. For the audio signal application example in the previous section, it may be decided that signal levels below 40dB will cause insignicant aliasing. The anti-aliasing lter used may have a bandwidth of 20kHz but the response is 40dB down starting from 24kHz. This means that the minimum sampling frequency has to be increased to 48kHz instead of 40kHz for the ideal lter. Alternatively, if we x the sampling rate, then we need an anti-alias lter with a sharper cut-off. Using the same audio example, if we want to keep the sampling rate at 44.1kHz, the anti-aliasing lter will need to have an attenuation of 40dB at about 22kHz. With a bandwidth of 20kHz, the lter will need a transition from 3dB at down to 40dB within 2kHz. This typically means that a higher order lter will be required. A higher order lter also implies that more components are needed for its implementation.
108
5.1.4 Practical Limits on Sampling Rates

As discussed in previous sections, the practical choice of sampling rate is determined by two factors for a certain type of input signal. On one hand, the sampling theorem imposes a lower bound on the allowed values of the sampling frequency. On the other hand, the economics of the hardware imposes an upper bound. This economics includes the cost of the analog-to-digital converter (ADC) and the cost of implementing the analog anti-alias lter. A higher speed ADC will allow a higher sampling frequency but may cost substantially more. However, a lower sampling frequency will put a more stringent requirement on the cut-off of the anti-aliasing lter, necessitating a higher order lter and a more complex circuit which again may cost more. In real-time applications, each sample is acquired (sampled), quantized and processed by a DSP. The output samples may need to be converted back to analog form. A higher sampling rate will means that there are more samples to be processed within a certain amount of time. If Tproc represents the total DSP processing time, then the time interval between samples Ts will need to be greater than Tproc . Otherwise the processor will not be able to keep up. This means that if we increase the sampling rate we will need a higher speed DSP chip.
5.2
Reconstruction of Analog Signals
A lot of applications requires an analog signal output. If the processing is done digitally, then the processed discrete-time signal will have to be converted to analog form. Since the information about the analog signal is contained in one of the replicas of the spectrum of the discrete-time signal, we need to isolate only one of those replicas. The easiest way to do so is to lter the discrete-time signal with a low-pass lter that has a passband in the Nyquist interval. There is an alternative way the Uniform Sampling Theorem can be stated: If the highest frequency contained in an analog signal x (t) is W and the signal is sampled uniformly at a rate fs of 2W or higher, then x (t) can be exactly recovered from its sample values using the formula
x (t) =
n=
x (nT ) hr (t nT )
(5.17)
Reconstruction of Analog Signals where hr (t) = sin (fs t) fs t 1 T = fs
109
(5.18) (5.19)
This is referred to as ideal reconstruction.
5.2.1
Ideal Reconstructor
The ideal reconstructor is simply an ideal low-pass lter with a cutoff frequency fs /2 where fs is the sampling frequency of the discrete-time signal. Its frequency response is given by T, |f | fs /2 Hr (f ) = (5.20) 0, otherwise where T = 1/fs . The impulse response is therefore the sinc function hr (t) = sinc Example 5.3. An analog signal x (t) = cos (10t) + cos (50t) + cos (90t) where t is in milliseconds, is preltered by an analog lter H (f ), producing y (t). The output y (t) is sampled at 40kHz, giving y (nT ). This discrete-time signal is passed through an ideal reconstructor that produces ya (t). Determine ya (t) when (i) H (f ) = 1; (ii) H (f ) is an ideal prelter with cutoff frequency 20kHz; and (iii) H (f ) is a practical prelter with a at magnitude response from DC to 20kHz and then attenuates monotonically at a rate of 60 dB/octave. Solution: (i) H (f ) = 1 implies that there is no preltering. Now, f1 = 5kHz f2 = 25kHz f3 = 45kHz t T (5.21)
110
Sampling at fs = 40kHz will cause aliasing for f2 and f3 and the aliased frequencies are f2a = 25 40 = 15kHz f3a = 45 40 = 5kHz So the reconstructed signal is ya (t) = cos (10t) + cos (30t) + cos (10t) = 2 cos (10t) + cos (30t) (ii) The prelter cuts off all frequency components above 20 kHz. Therefore, y (t) = cos (10t) and it is reconstructed without error with ya (t) = cos (10t) (iii) Recall that 1 octave is a doubling in the frequency. So the response at 40 kHz is 60 dB below that at 20 kHz. We have log2 log2 25 20 45 20 = 0.322 = 1.170
Therefore, attenuation at 25 kHz and 45 kHz are A2 = 60 (0.322) = 19.3dB A3 = 60 (1.170) = 70.1dB Amplitudes at these frequencies become |H (f2 )| = 10A2 /20 = 0.1084 |H (f3 )| = 10A3 /20 = 3.126 104 So ya (t) = cos (10t) + 0.1084 cos (30t) + 3.126 104 cos (10t) = 1.0003126 cos (10t) + 0.1084 cos (30t)
111
5.2.2 Staircase Reconstructor

While the ideal reconstructor is useful for analysis purposes, it cannot be implemented in practice. In reality, we need to rst convert the sequence of numbers that is the digital signal into some analog form. This analog signal will then be ltered to isolated only the frequency components in the Nyquist interval. A practical reconstructor is the staircase reconstructor. It generates a staircase approximation in response to each digital sample. Its impulse response is hr (t) = It can be expressed as hr (t) = u (t) u (t T ) This gives us the frequency response Hr () = 1 ejT j j 1 = 1 ejT j ejT /2 jT /2 = e ejT /2 j sin (T /2) jT /2 = e (5.23) 1, 0 t T 0, otherwise (5.22)
(5.24)
and so
sin (f T ) jf T e (5.25) f T The magnitude response of the staircase reconstructor is shown in Figure 5.8 together with the ideal reconstructor. H (f ) = T The staircase reconstruction is performed by Digital-to-Analog Converters (DACs).
5.2.3
Anti-image Postlters
The signal at the output of the staircase reconstructor yr (t) has a spectrum that is the product of the spectrum of the discrete-time signal and the frequency response of the
112
Staircase Reconstructor
Ideal Reconstructor
Sampled Signal Spectrum 0.8
Magnitude
0.6
0.4
0.2
0 2.5
1.5
0.5 0 0.5 Normalized Frequency (f/fs)
1.5
2.5
Figure 5.8: Magnitude Response of the Staircase and the Ideal Reconstructors. reconstructor. From Figure 5.8 it is obvious that the magnitude spectrumm of yr (t) will consists of most of the central replica together with some attenuated portions of other replicas. In order to obtain a good bandlimited reconstructed signal, yr (t) should be lowpass ltered so that only the central replica is retained. This low-pass lter, with a cutoff frequency of fs /2, is called an anti-image postlter. If the sampling frequency is near critical, then we need to have an anti-image postlter that has a very sharp roll-off because the transition band is narrow. Alternatively, if the sampling rate is high, then the replicas are far apart and so the roll-off will be gradual. Sharp roll-off implies high order lters have to be used. Since these are analog lters, they are typically implemented using op-amps, resistors and capacitors. Higher order lters involve more of these components and are therefore more expensive and more difcult to maintain. Lower order lters are therefore preferred. However, the tradeoff for using a lower order lter is to have a higher sampling rate and so we need a faster DAC and also a faster processor that can produce more samples within a given time. Real

15 dB/ octave
113
f (kHz) 12 8 4 0 4 8 12
Figure 5.9: Magnitude Spectrum of Sampled Signal in Example 5.4. designs therefore need to consider these issues. Example 5.4. The spectrum of a certain analog signal is at up to 4 kHz. Beyond 4 kHz, the spectrum attenuates at a rate of 15dB per octave. This signal is sampled at 12 kHz. What would be the amount of aliasing if no preltering is used? Solution: The spectrum of the sampled signal looks like that shown in Figure 5.9. By symmetry, x = y. Since 8 kHz is exactly one octave from 4 kHz, x = y = 15dB. So the aliased components are 15 dB or below the signal level. Example 5.5. For the signal given in Example 5.4, we now wish to suppress the highest aliased components within the 4 kHz band by more than 50dB. Determine the stopband attenuation required for the prelter. Solution: Referring to Figure 5.9, attenuation at 8 kHz is 15 dB. Since this is the highest aliasing component within the signals bandwidth, we need to increase the attenuation to 50 dB. In other words, the prelter must introduce a further attenuation of (50 15) = 35dB. Therefore the speccations for the prelter is a at passband up to 4 kHz, stopband edge at 8 kHz with minimum stopband attenuation of 35 dB.
114
Example 5.6. The sampled signal in Example 5.5 is ltered by a digital low-pass lter (which can be considered ideal) with cutoff frequency of 2 kHz. The ltered digital signal is reconstructed by a staircase DAC followed by a postlter. The overall reconstructor is required to suppress spectral images by more than 80dB. What is the attenuation requirement for the postlter? Solution: Since the bandwidth of the sampled signal is reduced to 2 kHz. The edge of the passband of the rst replica starts at 10 kHz instead of 8 kHz. Total attenuation is the sum of that introduced by the staircase reconstructor and the postlter (in dB), A (f ) = ADAC (f ) + AP OST (f ) The attenuation provided by the reconstructor is ADAC (f ) = 20 log10 = 20 log10 |H (f )| |H (0)| sin (f /fs ) f /fs
from (5.25). At f = 10kHz, ADAC (f ) = 14.4dB. Therefore the attenuation required of the postlter is given by AP OST (f ) = 80 14.4 = 65.6dB starting from f = 10kHz.
5.3 Quantization
5.3.1 Sample-and-Hold
Apart from sampling (discretization in time), the process of converting an analog signal into digital form also involve the discretization of the sampled signal amplitude or quantization. In practice, because the quantization process takes a nite amount of time, the sampled signal amplitude has to be held constant during this time. The sampling process is usually performed by a sample-and-hold circuit which can be logically represented
Quantization
115
as in Figure . The quantization process is performed by the analog-to-digital converter (ADC). The hold capacitor hold the sampled measurement of the analog signal x(nT ) for at most T seconds during which time a quantized value xQ (nT ) is available at the output of the analog-to-digital converter, represented as a B-bit binary number. The sample-and-hold and the ADC may be separate modules or may be integrated on the same chip. Typically the very fast ADCs require an external sample-and-hold device.
5.3.2 Uniform Quantization

The ADC assumes that the input values cover a full-scale range, say R. Typical values of R are between 1 to 15 volts. Since the quantized sampled value xQ (nT ) is represented by B-bits, it can take on only one of 2B possible quantization levels. If the spacing between these levels is the same throughout the range R, then we have a uniform quantizer. The spacing between quantization levels is called quantization width or quantizer resolution. For uniform quantization, the resolution is given by Q= R 2B (5.26)
The number of bits required to achieve a required resolution of Q is therefore B = log2 R Q (5.27)
Most ADCs can take bipolar inputs which means the sampled values lie within the symmetric range R R x(nT ) < (5.28) 2 2 For unipolar inputs, 0 x(nT ) < R (5.29) In practice, the input signal x(t) must be preconditioned to lie within the full-scale range of the quantizer. Figure 5.10 shows the quantization levels of a 3-bit quantizer for bipolar inputs.
116

Output Amplitude
A3 A2 A1 R/2 A 1 A 2 A 3 A 4 R/2
Input Amplitude
Figure 5.10: A Uniform 3-bit Quantizer Transfer Function.
5.3.3
Quantization Noise
Quantization error (or noise) is the difference between the actual sampled value and the quantized value. Mathematically, this is e(nT ) = x(nT ) xQ (nT ) or equivalently, e [n] = x [n] xQ [n] (5.31) Based on these equations, the quantization noise is modelled as an additive noise that is uncorrelated with the input value. This is not quite true unless the quantizer resolution is high since obviously the noise depends on the value of the input data as shown in the equations above. However, this assumption is usually made to simplify the analysis. Quantization noise is also normally assumed to be uniformly distributed within the quantization width. If x [n] lies between two quantization levels, it will either be rounded up or truncated. Rounding replaces x [n] by the value of the nearest quantization level. Truncation replaces x [n] by the value of the level below it. (5.30)
Quantization For rounding, the error is given by whereas for truncation, it is 0e<Q Q Q <e< 2 2
117
(5.32)
(5.33)
It is obvious that rounding produces a less biased representation of the analog values. The average error is given by
Q/2
e=
Q/2
ep(e)de
(5.34)
where p(e) is the distribution of e[n]. Assuming uniform distribution for quantization noise, i.e. p (e) = 1/Q, we have e= 1 Q
Q/2
ede = 0
Q/2
(5.35)
This means that on average half the values are rounded up and half down.
5.3.4
Signal-to-Quantization Noise Ratio
The mean-square value of the error gives us an idea of the average power of the error signal. It is given by
Q/2 2 e =
e2 p (e) de 1 Q
Q/2 Q/2
(5.36)
e2 de
Q/2 Q/2 Q/2
1 e3 = Q 3 Q2 = 12
(5.37)
118
2 For signal power x , the signal-to-quantization-noise ratio is therefore
SQN R = 10 log10 = 10 log10 = 10 log10
2 x 2 e 2 x Q2 /12 2 x (12) 2B 2 R
(5.38)
(5.39) R x (5.40)
= 20B log10 2 + 10 log10 12 20 log10 = 6.02B + 10.79 20 log10 R x
by using (5.37) and (5.26). Thus if we increase the number of bits of the ADC by one, the signal to quantization noise ratio improves by 6 dB. The above equation (excluding the bias value 10.79) gives us the dynamic range of the quantizer. Example 5.7. The dynamic range of the human ear is about 100 dB. If a digital audio system is to match this dynamic range, how many bits would be required to represent each signal sample? Ans: A 16-bit quantizer will provide a dynamic range of 96 dB.
5.4 Sampling Rate and Quantization Tradeoff

5.4.1 Power Spectral Density
The power spectral density (PSD) of a signal (or noise) is dened as the power per unit frequency. Thus the PSD of quantization noise is See =
2 e fs
(5.41)
since the spectra of sampled signals are unique only within a range equal to fs . The noise power within any particular range of frequencies [fa , fb ] can therefore be calculated by Pab = (fb fa ) See 2 (fb fa ) = e fs (5.42) (5.43)
Sampling Rate and Quantization Tradeoff
119
fs2 /2
fs1 /2
Figure 5.11: Power Spectral Density of Quantization Noise at Two Different Sampling Rates.
5.4.2 Oversampling
The total amount of quantization noise is dependent on the resolution of the quantizer. Given a quantizer, the total quantization noise power is the same. Since the power spectral density is a function of the sampling frequency fs , sampling at a higher rate will spread the total quantization noise power over a wider frequency range. Figure illustrates this fact. The total area under each distribution is the same. Let the quantization noise PSD be See1 when the sampling rate is fs1 , giving a total noise 2 power of e1 . At a higher sampling rate fs2 , the PSD becomes See2 for a total noise power 2 e2 . Since the total noise power is the same for any given quantizer,
2 2 e1 = e2
Thus See1 fs1 = See2 fs2 See1 fs2 = =L See2 fs1 where L is called the oversampling ratio. or (5.45) (5.46)
Now, suppose we have two quantizers, one is B1 bits and the other is B2 bits, with B2 < B1 . We desire to keep the output of these two quantizers to be of the same quality.
f s1 /2
e1 /f s1
e2 /f s2 f f s2 /2
(5.44)
120
In other words, the quantization noise level should be the same. One way to achieve this is to match each quantizer with different sampling rates: fs1 and fs2 (> fs1 ). In this case, we are keeping the PSD constant, i.e. See1 = See2 . So
2 2 e1 = e2 fs1 fs2 2 2 e1 = e2
(5.47) fs1 fs2 (5.48) (5.49)
= Now,
2 e1 = 2 e2
2 e2 L
Q2 1 12 Q2 2 = 12
(5.50) (5.51)
and Q1 = R2B1 Q2 = R2B2 So (5.49) becomes L =

2 e2 2 e1 Q2 2 = = 22(B1 B2 ) = 22B Q2 1
(5.52) (5.53)
(5.54) (5.55)
or B = 0.5 log2 L (5.56) Thus, by doubling the sampling rate (L = 2), we can reduce half a bit for the quantizer while maintaining the same quality of output. (5.56) is a tradeoff between sampling rate and quantization resolution. Example 5.8. CD quality sound requires 16-bit quantization. If we decided to use a 1-bit quantizer, what should be oversampling ratio be? Ans: L = 230 .

Highpass Filtered Noise Spectrum
121
noise remaining in signal band
e2 /f s2
fs2 /2
fs1 /2
f s1 /2
f s2 /2
Figure 5.12: Power Spectrum of Noise Before and After Noise Shaping.
5.4.3
Noise Shaping
The tradeoff of half a bit for doubling the sampling rate is too uneconomical. We can improve on that by using noise shaping techniques. Note that the quantization noise is white, i.e. its power spectrum is at. If the signal spectrum occupies just a portion of fs , then we can shape the noise spectrum so that it is lower within the signal band and higher for the rest of the frequencies by using a high-pass lter. The effect of the noise-shaping lter on the noise spectrum is shown in Figure . We can ensure that the signal bandwidth is small compared with the sampling frequency by oversampling. Consider the frequency response of a p-th order noise shaping high-pass lter HN S (f ) for |f | 2f fs
2p
(5.57)
|fs /2|. The average noise power within |f | |fs1 /2| is given by
2 e = 2 e2 fs2 fs1 /2 fs1 /2
2f fs
fs1 /2
2p
df
(5.58) (5.59) (5.60)
2 = e2 (2)2p 2 = e2
f 2p 2p+1 df fs1 /2 fs2 2p 1 2p+1 2p + 1 L
122 Therefore, B = p+ 1 2
log2 L
1 log2 2
2p 2p + 1
(5.61)
By using oversampling and noise shaping, doubling the sampling frequency gives us a 1 saving of p + 2 bits. This concept is being used in sigma-delta ADCs.
5.4.4
Non-uniform Quantization
One of the assumptions we have made in analysing the quantization error is that the sampled signal amplitude is uniformly distributed over the full-scale range. This assumption may not hold for certain applications. For instance, speech signals are known to have a wide dynamic range. Voiced speech (e.g. vowel sounds) may have amplitudes that span the whole full-scale range while softer unvoiced speech (e.g. consonants such as fricatives) usually have much smaller amplitudes. Also, an average person only speaks 60% of the time while she/he is talking. The remaining 40% is silence with a negligible signal amplitude. If uniform quantization is used, the louder voiced sounds will be adequately represented. However, the softer sounds will probably occupy only a small number of quantization levels with similar binary values. This means that we would not be able to distinguish between the softer sounds. As a result, the reconstructed analog speech from these digital samples will not nearly be as intelligible as the original. To get around this problem, non-uniform quantization can be used. More quantization levels are assigned to the lower amplitudes while the higher amplitudes will have less number of levels. This quantization scheme is shown in Figure 2-17. Alternatively, a uniform quantizer can still be used, but the input signal is rst compressed by a system with an input-output relationship (or transfer function) similar to that shown in Figure 5.14. The higher amplitudes of the input signal are compressed, effectively reducing the number of levels assigned to it. The lower amplitude signals are expanded (or non-uniformly amplied), effectively making it occupy a large number of quantization levels. After processing, an inverse operation is applied to the output signal (expanding it). The system that expands the signal has an input-output relationship that is the inverse of the compressor. The expander expands the high amplitudes and compresses the low amplitudes. The whole process is called companding (COMpressing and exPANDING).

Output
123
Input
Figure 5.13: Non-uniform Quantization. 5.4.4.1 Companding Companding is widely used in public telephone systems. There are two distinct companding schemes. In Europe, A-law companding is used and in the United States, -law companding is used. -law compression characteristic is given by the formula: ln 1 + y = ymax where sgn (x) =
|x| xmax
ln (1 + )
sgn (x)
(5.62)
+1, x 0 1, x < 0
(5.63)
Here, x and y represent the input and output values, and xmax and ymax are the maximum positive excursions of the input and output, respectively. is a positive constant. The North American standard species to be 255. Notice that = 0 corresponds to a linear input-output relationship (i.e. uniform quantization). The compression characteristic is shown in Figure 5.14.
124
The A-law compression characteristic is given by |x| y A( xmax ) sgn (x) , 1 0 < x|x| A max 1+ln A max y= |x| y 1+ln[A( xmax )] sgn (x) , 1 < |x| < 1 max 1+ln A A xmax
(5.64)
Here, A is a postive constant. The European standard species A to be 87.6. Figure 5.15 shows the characteristic graphically.
1 0.9 u=255 0.8 0.7 Compressed Input 0.6 0.5 0.4 0.3 0.2 0.1 0 u=10 mulaw compression
u=100
0.1
0.2
0.3
0.4
0.5 0.6 Normalized Input
0.7
0.8
0.9
Figure 5.14: -law Compression Function
5.5
Dithering
Another assumption we have made in analysing quantization noise is that it is assumed to be uniformly distributed over the quantization width. If the noise is not uniformly distributed, quantization distortion results.
Dithering
125
Figure 5.15: A-law Compression Function We shall illustrate quantization distortion through an example. A low amplitude sinusoid is being sampled and quantized. The samples of the sinusoid are given by x [n] = A cos (2f0 n) where A is less than the quantization resolution. Let fs = 40 samples per second and A = 0.75Q (5.66) So for a 1 kHz sinusoid, the actual sampling rate is 40 kHz. Figure 5.16 shows the original and the quantized signals. Note that the quantized signal only occupies three of the available quantization levels. The frequency spectrum of this quantized signal is shown in Figure 5.17. It has peaks at f0 , and the odd harmonic frequencies 3f0 , 5f0 , etc. Clearly the odd harmonics are artefacts of the quantization process and can be considered as the spectrum of the quantization noise signal which, in this case, is not white. (5.65)
126

Undithered Quantization
1 0.8 0.6 0.4 0.2 Units of Q 0 0.2 0.4 0.6 0.8 1
100
200
300
400
500 Index n
600
700
800
900
1000
Figure 5.16: Original and Quantized Signals This problem can be overcome by adding a dither v(n) to the original sampled signal so that y(n) = x(n) + v(n). Various types of dither can be added. Two of them which are of practical interest are rectangular and triangular dither. They are so called because the distribution of the random signal samples are rectangular and triangular in shape respectively. The distributions are shown in Figure 5.18. The addition of dither to the original signal will increase its average quantization noise power. Recall that the average noise power for uniform quantization is Q2 /12. The addition of rectangular dither will double this average noise power and the addition of triangular dither will triple it. However, if we look at the frequency spectrum of the dithered and quantized signal of the example we have been considering (Figure ), we will notice that the noise spectrum now appears to be white and the odd harmonic artefacts are not there any more.
Summary
127
Figure 5.17: Spectrum of the Quantized Signal It must be emphasized that in general, the sampling process will cause all the odd harmonics that lie outside of the Nyquist interval (out-of-band harmonics) to be aliased back into the interval (in-band non-harmonic frequencies). So the overall spectrum will contain peaks at frequencies other than the odd harmonics.
5.6
Summary
While the uniform sampling theorem is the basis for choosing sampling rates, we have shown that in practice, issues such as the complexity of the anti-aliasing prelter and the anti-image postlter will need to be considered. The sampling theorem merely provide us with a lower bound on the sampling rate. The upper bound is limited by the speed at which the processor can handle the samples and the complexity of the processing required.
128
p(v)
p(v)
1/Q
1/Q
Q/2
Q/2
(a) Uniform distribution
(b) Triangular Distribution
Figure 5.18: Amplitude Distributions of Rectangular and Triangular Dither
Figure 5.19: The Dithered Signal and Its Quantized Version
Summary
129
Figure 5.20: Spectrum of Quantized Signal with Dithering One can reconstruct almost perfectly an analog signal that has been sampled fast enough. However, once the sampled signal is quantized, it is forever distorted. The amount of distortion can be minimized by increasing the resolution of the quantizer. Increasing quantizer resolution implies increasing the wordlength of the data to be processed. There is a tradeoff between the sampling rate and the quantizer resolution for a certain the quality of the output. All these are practice design issues to be addressed for each DSP system that involves analog-to-digital and digital-to-analog conversions.
130
Chapter 6 Digital Filters and Their Realizations

In this chapter, we are concern about how the output of a digital lter is computed given its coefcients and the input signal. Even though the digital lter equation is simple, there are many different ways by which the output could be computed and each method has its own characteristics. We shall also study the consequences of representing the coefcients and the results of the computing with a nite number of bits. But rstly, a brief review of Z transform is given. We assume that the readers are already somewhat familiar with Z transform. Those who are not should consult one of the relevant references.
6.1
A Review of Z Transform
6.1.1 Denition
The Z transform of a discrete-time sequence x[n] is dened by
X (z) =
n=
x[n]z n
(6.1)
if the series converges. Here, z is a complex variable. The region of convergence (ROC) of a Z-transform is the range of z for which the transform converges. Since z is complex, the ROC is a region on the complex plane. X (z) and its ROC uniquely determines the sequence x[n].
132 Signal x [n] [n] u [n] an u [n] nan u [n] cos (0 n) u [n] an u [n 1] Z Transform 1
1 1z 1 1 1az 1 az 1 (1az 1 )2 1z 1 cos 0 12z 1 cos 0 +z 2 1 1az 1
Digital Filters and Their Realizations Region of Convergence All z |z| > 1 |z| > |a| |z| > 1 |z| > 1 |z| < |a|
Table 6.1: Some Z-Transform Pairs Example 6.1. Determine the Z-transform of the sequence x [n] = (0.5)n u [n] where u [n] is the unit step sequence. Example 6.2. Determine the Z-transform of the sequence x [n] = (0.5)n u [n 1] Some common Z-transform pairs and their ROCs are listed in Table 6.1.
6.1.2
Properties
The Z-transform is linear and so the superposition principle holds. a1 x1 [n] + a2 x2 [n] a1 X1 (z) + a2 X2 (z) (6.2)
Linearity
where a1 and a2 are arbitrary constants and x1 [n] X1 (z) x2 [n] X2 (z) Delay If x [n] X (z) then a delay of D samples x [n D] z D X (z)
A Review of Z Transform
133
Convolution Convolution in the time-domain will be transformed to multiplication in the Z domain. y [n] = x [n] h [n] Y (z) = H (z) X (z) (6.3) Causality Consider a causal sequence x [n] = pn u [n] + pn u [n] + 1 2 with Z transform X (z) = The ROC is given by |z| > max |pi |
i
(6.4)
1 1 + + 1 1 p1 z 1 p2 z 1
(6.5) (6.6)
Causal signals are characterized by ROCs that are outside the maximum pole circle. If a sequence x[n] is causal, then only negative powers of z appear in X (z). If x[n] is stricly anti-causal, then X (z) only consists of positive powers of z. Stability The necessary and sufcient condition for a sequence x [n] to be stable is that its Z-transform X (z) has a ROC which contains the unit circle. For causal signals, (6.6) holds. That means a causal and stable sequence must have all its poles inside the unit circle. Similarly, a causal and stable system must have a transfer function with all the poles inside the unit circle.
6.1.3
Inverse Z Transform
The inverse Z-Transform is dened by x [n] = 1 2j X (z) z n1 dz

C
(6.7)
with the contour integral evaluated in the anti-clockwise direction. Example 6.3. Find the inverse Z-transform of X (z) = 2 2.05z 1 1 2.05z 1 + z 2
134
Digital Filters and Their Realizations
6.2
System Transfer Functions
6.2.1 Poles and Zeros

Consider a rational Z-transform with real-valued coefcients ak and bk : N (z) b0 + b1 z 1 + + bM z M H (z) = = D (z) a0 + a1 z 1 + + aN z N which can be factorized as H (z) = b0 N M (z z1 ) (z z2 ) (z zM ) z a0 (z p1 ) (z p2 ) (z pN ) (6.9) (6.8)
H (z) is zero for z = z1 , z2 , . . . , zM . These are the zeros of H (z). The transfer function is innite at z = p1 , p2 , . . . , pN . These are the poles of H (z). The poles and zeros may be real-valued or complex-valued. In the case where they are complex, they must occur in complex conjugate pairs since ak and bk are real-valued. If a0 = 1 and all other ak are zeros, then H (z) is an all-zero system. We say that it has M zeros and an M -th order pole at z = 0. The transfer function now becomes
M
H (z) =
k=0
bk z k
(6.10)
The impulse response sequence is therefore given by {b0 , b1 , . . . , bM } and is nite if M is nite.
6.2.2
Frequency Response
The frequency response of an LTI system with a rational transfer function H (z) can be obtained by H () = H (z) |z=ej = =
M jk k=0 bk e N jk k=0 ak e M j ) b0 k=1 (1 ck e N j ) a0 k=1 (1 dk e
(6.11) (6.12) (6.13)
System Transfer Functions The gain is the log magnitude of H (). It is usually expressed in decibels (dB): Gain (dB) = 20 log10 |H ()| The phase response is dened as () = H ()
135
(6.14)
(6.15)
Group Delay Consider a system with impulse response h [n] = [n nd ] which consists of only one impulse at some arbitrary index nd . The frequency response of this system is given by
H () = = e
h [n] ejn
n= jnd
(6.16) (6.17)
for || . The phase response is therefore H () = () = nd We dene the phase delay as P hase Delay = () (6.19)
(6.18)
In this case it is equal to nd which is independent of . We call systems with phase delays independent of linear phase systems. An alternative way to measure the linearity of the phase response is to use the group delay dened by d () g () = (6.20) d If H () has a linear phase response, then its group delay is a constant.
136
Minimum Phase Systems The stability and causality of a system are solely determined by the poles of the transfer function; they do not concern the zeros. For certain problems, we need to ensure that the inverse system is also causal and stable. In these cases, both zeros and poles of the transfer function must be inside the unit circle. This kind of systems are called minimum phase systems. Linear Phase FIR Systems FIR transfer functions can have exact linear phase. Since such systems have a constant group delay g , they must have a frequency response of the form H () = G () ejg +j (6.21) where G () is a real-valued function and g and are constants. The impulse response of this lter is given by the inverse DTFT: h [n] = = = =
1 H () ejn d 2 1 G () ejg +j ejn d 2 ej G () ej(ng ) d 2 ej g [n g ]
(6.22)
where g [n] is the IDTFT of G (). So h [n + g ] = ej g [n] Since G () is real, g [n] = g [n] and so from (6.23) we have ej h [n + ] = ej h [n + g ] or h [n + g ] = ej2 h [n + g ] (6.25) If h [n] is causal and of nite length (0 n N 1), then the group delay must be g = (N 1) 2 (6.26) (6.24) (6.23)
Digital Filter Realization and so (6.25) becomes h [n] = ej2 h [N n 1]
137
(6.27)
for 0 n N . For most cases, we design lters with real-valued h [n] and therefore ej2 must be real, implying k = 2 for integer values of k. So (6.27) can be rewritten as h [n] = (1)k h [N n 1] (6.28)
This means that the lter impulse response must be either symmetric or anti-symmetric. In general, we can identify four different types of linear phase responses: k is even and N 1 is even (or N is odd). k is even and N 1 is odd (or N is even). k is odd and N 1 is even (or N is odd). k is odd and N 1 is odd (or N is even).
6.3 Digital Filter Realization

Given the transfer function of a digital lter, we have all that we needed to compute the output of the lter given some input. However, even though the equation looks simple, there are many ways by which the output could be computed. When we talk about realization, we are talking about a structure that relates the input and output of the lter. The most common way to specify a particular realization, or lter structure, is through the use of block diagrams. The blocks within these block diagrams can be implemented either as a piece of digital hardware or as a program to be executed by a DSP chip. We shall describe several lter structures. Considerations for the choice of a lter structure include ease of programming on a particular DSP chip, and the regularity of the VLSI (Very Large Scale Integrated) design. Some structures are more sensitive to (quantization) errors in the coefcients. In some cases, such as IIR lters, the stability of the lter may depend on an appropriate realization.
138
b0 x[n]
+
z 1 b1 a1 z 1 y[n1] z 1 b2 a2 y[n2] z 1
y[n]
x[n1]
x[n2]
Figure 6.1: Direct Form Realization of a Second Order IIR Filter.
6.3.1
Direct Form I
Consider a simple second order IIR lter with transfer function H (z) = B (z) b0 + b1 z 1 + b2 z 2 = A (z) 1 + a1 z 1 + a2 z 2 (6.29)
The input and output samples are related by y [n] = a1 y [n 1] a2 y [n 2] + b0 x [n] + b1 x [n 1] + b2 x [n 2] (6.30)
Direct form realization is simply a realization based on the direct implementation of this difference equation. This is illustrated in Figure 6.1. In the gure, the z 1 symbol represents the delay of one sample. In actual implementations, it would represent shift registers or a memory location in RAM. The three basic elements of the structure are illustrated in Figure 6.2. This structure involves 4 unit delays, 5 multiplications and 4 additions. The structure also requires the storage of past values of x and y. In this case, 4 memory units or registers, excluding those for the immediate input and output, are required. In general, for an IIR transfer function M k k=0 bk z H (z) = (6.31) 1 + N ak z k k=1
Digital Filter Realization
139
+ ....
(a) Adder
z1
(b) Multiplier
(c) Unit Delay
Figure 6.2: Three basic elements in the realization. we require (M + N + 1) multiplications, (M + N ) additions and (M + N ) memory units. Notice that the right hand side of (6.30) consists of two main operations: multiplications and additions. Each sample of the output y or input x is multiplied by a coefcient. The result is then stored or accumulated for addition. We call these two basic operations of the direct form structure Multiply-and-Accumulate (MAC). Dedicated DSP processors have special instructions for MAC that are very efcient. Figure 6.1 also illustrates the two parts of the lter structure. All the numerator terms shown on the left hand side of the adder block are the feed-forward elements. The denominator terms that depend on the previous output samples are feeding back. This realization is called direct form I or simply direct form.
6.3.2
Canonical Form or Direct Form II
We can view the structure in Figure 6.1 as a cascade of two sub-systems H1 (z)and H2 (z) as shown in Figure 6.3 where H1 (z) = B (z) 1 H2 (z) = A (z) (6.32) (6.33)
Note that H1 (z) consists of only feed-forward terms and H2 (z) only has feedback terms. Since this is a linear system, reversing the order of the cascade will not change the overall
140
H (z)
1

H2 (z) 1/A(z)
B(z)
b0 x[n]
+
z 1 b1
+
a1 z 1 y[n1]
y[n]
x[n1]
z 1 b2 x[n2] a2
z 1 y[n2]
Figure 6.3: Direct Form I Structure as a Cascade of Two Sub-systems. system transfer function H (z). Figure 6.4 depicts the cascade of the sub-systems with the order interchanged. The output of the sub-system H2 (z) is now the input to the subsystem H1 (z) . The two sets of delays in the middle operate on the same data w[n]. Therefore we do not need two separate sets of delays; they can be merged into one as shown in Figure 6.5. This realization is known as the Direct Form II realization or the Canonical Form realization. It is not difcult to see that the canonical form implements the original IIR lter equation (6.30). We have, from the adder on the left hand side of Figure 6.1,
2
w[n] =
k=1
ak w [n k] + x [n] (6.34)
= a1 w [n 1] a2 w [n 2] + x [n]

H (z)
1
141
H2 (z) 1/A(z)
B(z)
x[n]
+
z 1 a1
w[n]
b0
y[n]
z 1 b1
z 1 a2
z 1 b2
Figure 6.4: Interchanging the Order of the Sub-systems.
The output of the adder on the right is given by
y [n] =
k=0
bk w [n k] (6.35)
= b0 w [n] + b1 w [n 1] + b2 w [n 2]
Substituting the expressions for w [n], w [n 1] and w [n 2] that can be obtained from
142 (6.34) into (6.35), we have
y [n] = b0 (a1 w [n 1] a2 w [n 2]) + b0 x [n] +b1 (a1 w [n 2] a2 w [n 3]) + b1 x [n 1] +b2 (a1 w [n 3] a2 w [n 4]) + b2 x [n 2]
2
(6.36)
=
k=0
bk x [n k] a1 (b0 w [n 1] + b1 w [n 2] + b2 w [n 3]) a2 (b0 w [n 2] + b1 w [n 3] + b3 w [n 4])

2
(6.37) (6.38)
=
k=0
bk x [n k] a1 y [n 1] a2 y [n 2]
which is the same as (6.30). The structure in Figure 6.5 requires 5 multiplications and 4 additions similar to direct form I. However, the number of memory units required is reduced to 2 (for w [n 1] and w [n 2]). In fact, the word canonical refers to the fact that this structure requires the least amount of storage. Referring to (6.31), the canonical form requires (M + N + 1) multiplications, (M + N ) additions, and max (M, N ) memory units, excluding the current input and output samples. Although we have been using the second order IIR lter as example, both direct form structures can easily be generalized to lter transfer functions of higher orders. Example 6.4. Given the digital lter transfer function 2 3z 1 + 4z 3 H (z) = 1 + 0.2z 1 0.3z 2 + 0.5z 4 derive the direct form I and II realizations. Relative Merits The canonical form is commonly used because of the following properties: 1. requires a minimum of storage space; and 2. good round-off noise property.
143
w[n] b0
x[n]
+
a1
y[n]
z 1 b1
z 1 a2 b2
Figure 6.5: Direct Form II or the Canonical Form. The disadvantage is that it is susceptible to internal numeric overow. The input to the lter can be scaled to avoid overows. On the other hand, the direct form I structure does not require scaling because there is only one adder. Thus if scaling is not desirable, direct form I realization is preferred.
6.3.3
Transposed Structures
The direct form I and II structures can be transposed. The transposed structures can be obtained by 1. reversing all the signal ow directions; 2. change node (connection points) into adders and adders into nodes; and 3. exchanging the input with the output. The transposed direct form I and II are shown in Figures 6.6 and 6.7 respectively. In both cases, the transfer function remains the same after transposition. Note that for the tranposed structures, the inputs x [n] and outputs y [n] are no longer stored. Only some computed intermediate values are stored. This allows for some optimization of the computation since one sample of x [n] is now multiplied with different coefcients.
144
y[n]
+
z
1
b0
+
z
1
x[n]
b1
a1
+
z
1
+
z
1
b2
a2
Figure 6.6: Transposed Direct Form I Realization of a Second Order IIR Filter.
6.3.4
FIR Filters
For an FIR lter, the denominator polynomial is simply equal to 1, i.e. A (z) = 1. So only the feed forward elements exist. The direct form realization is usually drawn in a slightly different way as shown in Figure 6.8. This corresponds to the FIR equation we have been using for a second order lter
2
y [n] =
k=0
h [k] x [n k]
(6.39)
where h [n] is the impulse response or lter coefcients. For linear phase FIR lters, the lter coefcients are symmetric or anti-symmetric. So for an N -th order lter, the number of multiplications can be reduced from N to N/2 for N even and to (N + 1)/2 for N odd. Figure 6.9 shows the direct form realization of an odd order linear phase FIR lter that takes advantage of this saving.

b0 x[n]
145
y[n]
+
z
a1
1
b1
+
z
1
a2
b2
Figure 6.7: Transposed Direct Form II Realization of a Second Order IIR Filter.
x[n]
x[n1]
x[n2]
h[0]
h[1]
h[2]
+
y[n]
Figure 6.8: Direct Form I Realization of an FIR Filter.
146
x[n] x[n1]
+
z
1
+
z
1
+
z
1
h[0]
h[1]
h[N/21]
h[(N1)/2]
y[n]
Figure 6.9: Direct Form Realization of Odd Order Linear Phase FIR Filter. It is interesting to note that a transposed FIR direct structure can also be obtained using the method discussed for transposing IIR structures. The resulting structure for a second order lter is shown in Figure 6.10.
6.3.5
Cascade Realization
A general transfer function H (z) of order N > 2 can be factorized into K factors so that
K
H (z) =
i=1 K
Hi (z) bi0 + bi1 z 1 + bi2 z 2 1 + ai1 z 1 + ai2 z 2
(6.40) (6.41)
=
i=1
where K = N/2 for even values of N and K = (N + 1) /2 if N is odd. We may view each second order section with transfer function Hi (z) as a subsystem of the whole
147
x[n]
h[0]
h[1]
h[2]
y[n]
Figure 6.10: Transposed FIR Filter Structure Realization of an FIR Filter.
x[n]
H0 (z)
H 1 (z)
HK1(z)
y[n]
H(z)=H 0 (z)H1 (z)...HK1(z)
Figure 6.11: Cascade of Transfer Functions
148
system H (z) and so the full system is made up of a cascade of these subsystems as depicted in Figure 6.11. Note that when we say second order section we really mean up to second order section. Some of the bi2 and ai2 can be zero. So the numerator and denominator polynomial orders can be less than $2K$. Also, the coefcients of each section are real-valued. Each second order section can be implemented using direct, canonical or transposed forms. Note that the second order transfer functions Hi (z) formed are not unique. Factors in the numerator and denominator polynomials can be paired in an arbitrary way to form Hi (z). But the overall transfer function H (z) remains the same. In practice, the pairing and ordering of the second order sections may affect the numeric accuracy of the resulting lter. The internal multiplication in each section may generate a certain amount of roundoff error. This error is then propagated to the next section. The round-off error of the overall output is different for each combination of second order sections. Naturally we want to achieve a minimum amount of round-off error. This is a difcult problem to solve. In practice, some trial-and-error would be needed. Fortunately, most IIR lters do not have a high order so the number of possible combinations are not too large. The following rules of thumb can be used for cascade realizations: Pair the numerator and denominator quadratic pairs with roots that are closest to one another. In this way, the poles and zeros tend to cancel one another and so the sections will not be producing output values that are very large or very small. Put the section with the denominator root having magnitudes that are closest to one as the last section. Systems with poles near the unit circle are the least stable ones. Example 6.5. Determine the cascade realization of the following lter: H (z) = Solution: Factorize H (z) into H (z) = (1 + z 1 ) (1 + z 1 ) (1 + z 1 ) (1 + z 1 ) = (1 0.5z 1 ) (1 0.25z 1 ) (1 0.25z 1 ) (1 0.5z 1 ) 1 + 2z 1 + z 2 1 0.75z 1 + 0.125z 2
The two sub-systems can be realized in any form.
149
6.3.6 Parallel Realization

An alternative to cascade form for higher order IIR lters is the parallel form. Here the transfer function H (z) is expanded as a sum of rst and second order terms.
N
H (z) = C +
k=1
Hk (z)
(6.42)
where C = Hk (z) = bN aN Ak 1 pk z 1
for an N -th order transfer function. In this case, the individual subsystem transfer functions are summed to form the overall transfer function. Thus the subsystems are connected in parallel in contrast with the cascade form. This is shown in Figure 6.12. The whole IIR lter now consists of a parallel bank of rst order lters. Both Ak and pk can be complex-valued. If pk is complex, then its complex conjugate will also be appear since the coefcients of H (z) are real-valued. We can combine the pair of complex conjugate terms to avoid having to deal with complex numbers. So the transfer function of the combined subsystems are second order sections given by bk0 + bk1 z 1 1 + ak1 z 1 + ak2 z 2
Hk (z) =
The advantage of the parallel form compared with the cascade form for IIR lters is that the ordering of the subsystems are unimportant since they are in parallel. Scaling is easier as it can be carried out for each sub-system independently. Furthermore, round-off errors in each block are not propagated to the next block. Example 6.6. Realize the lter given in Example 6.5 in parallel form. Solution:
150
C
H 1 (z)
H2 (z) x[n]
y[n]
H N (z)
Figure 6.12: Parallel Form Realization
H (z) =
1 + 2z 1 + z 2 1 0.75z 1 + 0.125z 2 7 + 8z 1 = 8+ (1 0.5z 1 ) (1 0.25z 1 ) 18 25 = 8+ 1 1 0.5z 1 0.25z 1
The realization using canonical form for each sub-system is shown in Figure.
6.3.7
Circular Buffers
Sample-by-sample processing of an IIR lter implemented in the canonical form can be expressed in pseudo-code as follows.
151
x[n]
+
z
0.5
1
18
y[n]
+
z
0.25
1
25
Figure 6.13: Parallel Canonical Form Realization of the Transfer Function in Example 6.6.
152 For each input sample x do: w(0)=x; for i=1,2,...,M do: w(0)=w(0)-a(i)w(i); y=b(L)w(L); for i=L-1,L-2,...,0 do: w(i+1)=w(i); y=y+b(i)w(i); end;
The number of multiply-accumulate (MAC) operations required is to compute each output sample y is M (L + 1). If each MAC can be performed in one instruction cycle with time Tinstr , then the total time required for the MACs will be TM AC = M (L + 1) Tinstr The algorithm also requires L data shifts. Similarly, for an (N + 1)-th order FIR lter, the sample processing can be expressed as: (6.43)
For each input sample x do: w(0)=x; y=b(N)w(N); for i=N-1,...,0 do:
Digital Filter Realization w(i+1)=w(i); y=y+b(i)w(i); end;
153
Thus to compute each output sample, (N + 1) MAC operations and N data shifts are required. The processing time if only MAC operations are considered is TM AC = (N + 1) Tinstr (6.44)
However, if each data shift also requires one instruction cycle, then the total processing time for the IIR lter is Tproc = M (L + 1) Tinstr + LTinstr and for the FIR lter it is Tproc = (N + 1) Tinstr + N Tinstr (6.46) (6.45)
Obviously the time for shifting the coefcients is a signicant portion of the total processing time. Since shifting data is much easier to perform than MAC, it should not have taken so much time. One way to avoid coefcient shifting altogether is to use a circular buffer. It is the equivalent of a circularly linked list. The input data samples are stored in the circular buffer. A pointer points to the starting location of the buffer. In order to illustrate how this works, let us consider a length-4 impulse response which requires a length-4 circular buffer. If xn is the current input data sample, then to compute the output sample, we need to perform: yn = h0 xn + h1 xn1 + h2 xn2 + h3 xn3 (6.47) where hi , i = 0, 1, 2, 3 is the impulse response and xn1 , xn2 , xn3 are the 3 past samples immediately preceding xn . The input samples are placed in the circular buffer as shown in Figure 6.14. The pointer is pointing at the current input sample. The next input sample xn+1 replaces the oldest sample xn3 since it is not longer needed to compute yn+1 . The pointer is updated to point at the location of xn+1 . The output yn+1 is now computed by multiplying the data that the pointer is pointing to with h0 and the rest of the product terms in the same anti-clockwise direction.
154
h1 x n1

h2 Next sample instant x n1
h2
x n2
xn
h0
h3
x n2 p
xn
h1
x n3 h3
x n+1 h0
Figure 6.14: A Length 4 Circular Buffer. Notice that now we no longer need to shift the input data as we used to do with a linear array. Once a new data sample is placed in the right place in the circular buffer, an update of the pointer is all we need and the same procedures are executed. All digital signal processors now support circular buffer operations by having dedicated address generators. Once the size of the circular buffer and the addresses of the coefcients and input data are set, the operation is performed automatically. This give a signicant improvement in ltering performance, in addition to the single instruction MAC.
6.4
Finite Wordlength Effects
The number of bits that is used to represent numbers, called wordlength, is often dictated by the architecture of DSP processor. If specialized VLSI hardware is designed, then we have more control over the word-length. In both cases we need to tradeoff the accuracy with the computational complexity. No matter how many bits are available for number representation, it will never be enough for all situations and some rounding and truncation will still be required. A number within a digital system is always represented by a nite number of bits because
155
of the nite wordlength of data registers used. Quantizing a number to a xed number of bits has an inherent limitation on the accuracy of the representation. For FIR digital lters, the nite word-length affects the results in the following ways: 1. Coefcient quantization errors. The coefcient we arrived at in the approximation stage of lter design assume that we have innite precision. In practice, however, we have the same word-length limitations on the coefcients as that on the signal samples. The accuracy of lter coefcients will affect the frequency response of the implemented lter. 2. Round-off or truncation errors resulting from arithmetic operations. Arithmetic operations such as addition and multiplication often give result that require more bits than the word-length. Thus truncation or rounding of the result is needed. Some lter structures are more sensitive to these errors than others. 3. Arithmetic overow. This happens when some intermediate results exceed the range of numbers that can be represented by the given word-length. It can be avoided by careful design of the algorithm scale the results are appropriate stages. For IIR lters, our analysis of nite word-length effects will need to include one more factor: Product round-off errors. The round-off or truncation errors of the output sample at one time instant will affect the error of the next output sample because of the recursive nature of the IIR lters. Sometimes limit cycles can occur. We shall rst review how xed and oating point numbers are represented before going into the details of the above effects.
6.4.1
Fixed-point Representation
The general xed point format is basically the same as the usual familiar decimal representation of numbers. It consists of a string of digits with a decimal point. The digits to the left of the decimal point are the integer part and those to the right are the fractional part of the number. X = (bB bB1 b1 b0 b1 b2 bA+1 bA )r
B
(6.48) (6.49)
=
i=A
bi r i
156
for 0 bi r 1 where bi are the digits and r is the radix or base. Example 6.7. (12.34)10 = 1 101 + 2 100 + 3 101 + 4 102 (110.01)2 = 1 22 + 1 21 + 0 20 + 0 21 + 1 22 We shall focus on binary representations as this is the format we need to deal with in DSP. In this case, the digits are bits. So bA is the least signicant bit (LSB) and bB is the most signicant bit (MSB). Naturally the binary point (as opposed to decimal point) between b0 and b1 does not physically exist and so the user must know beforehand how to interpret a particular binary number. Non-negative integers can easily be represented by an n-bit pattern (B = n 1 and A = 0). Since we need to deal with fractional numbers, the fraction format (B = 0 and A = n) is normally used. This allows us to respresent numbers in the range 0 to 1 2n . This is because multiplication of two numbers that are less than 1 will give us a result that is less than 1. For signed numbers, an extra sign bit is required. So each number X is normalized to satisfy 1 X < 1 Positive fractions are given by X = 0.b1 b2 bn+1
n
(6.50) (6.51) (6.52)
=
i=1
bi 2i
Negative fractions can be represented in the twos complement form. That is, X = 1.b1 b2 bn+2 bn+1 1 where is the exclusive-OR operator. Example 6.8. 7 = 0.111 8 7 = 1.001 8 Multiplication of two xed point numbers each n bits in length will generally give a product which is 2n bits in length. The product therefore has to be truncated or rounded off to n bits, producing truncation or round-off errors. (6.53)
157
6.4.2 Floating-point Representation

Floating point representations are able to cover a much wider range of numbers for the same number of bits compared with xed point representations. They normally consist of a mantissa M , which is the fractional part of the number, and an exponent E, which can be either positive or negative. Hence a number X is given by X = M 2E (6.54)
with 1/2 M < 1. The mantissa and the exponent require their individual sign bits. Given a total number of bits available for representing a number, a number of different oating point formats can result. In the past, individual computer manufacturers used their own format for their own products. A common standard oating point format has been adopted by the Institute of Electrical and Electronic Engineers (IEEE) which is usually referred to as the IEEE 754 standard. It denes the way zero is represented, the choice of M and E, the handling of overows and other issues. For a 32-bit representation, the single precision oating point number is dened as X = (1)S 2E127 (M ) (6.55)
where S is the sign bit, E occupies 8 bits, and M is 23 bits long in a format shown in Figure 6.15. The following rules apply: 1. If E = 255 and M = 0, then X is not a number (denoted NaN). 2. If E = 255 and M = 0, then X is innity (denoted InF). 3. If 0 < E < 255, then X = (1)S 2E126 (1.M ). 4. If E = 0 and M = 0, then X = (1)S 2126 (0.M ). 5. If E = 0 and M = 0, then X is zero. Here (0.M ) is a fraction and (1.M ) is a number with one integer bit and 23 fractional bits.
158
Bits 0 1 S sign E 8 9

31 M mantissa
exponent
Figure 6.15: IEEE-754 Floating Point Format Example 6.9. What is the decimal value of the following IEEE-754 formatted oating point number? 0100000101010 00 Ans: 13. Floating point representations can naturally represent a much larger range of numbers than a xed point one with the same number of bits. However, it should be noted that the resolution does not remain the same throughout this range. This means that the distance between two oating point numbers increases as the number is increased. On the other hand, the resolution of xed point numbers is constant throughout the range. When two oating point numbers are multiplied, the mantissas are multiplied and the exponents are added. But if we want to add two oating point numbers, the exponents of the two numbers must be equal. The one with the small exponent is adjusted by increasing the exponent and reducing the mantissa. This adjustment could result in a loss in precision in the mantissa. Overow occurs in multiplication when the sum of the two exponents exceed the dynamic range of the representation for the exponent.
6.4.3
Coefcient Quantization
Let us consider the low-pass linear phase FIR lter with a normalized cutoff frequency of /4. The lter coefcients are quantized to 8 bits. The magnitude response of the lter, both before and after coefcient quantization, is shown in Figure 6.16. If the minimum stopband attenuation allowed is 48dB, then it is violated by the quantized lter. Clearly in this case, more than 8 bits are required for the lter coefcients.
159
The minimum number of bits needed for the lter coefcients can be found by computing the frequency response of the coefcient quantized lter. A trial-and-error approach can be used. However, it will be useful to have some guideline for estimating the word-length requirements of a specic lter. The quantized coefcients and the unquantized ones are related by hq [n] = h [n] + e [n] (6.56)
for n = 0, 1, . . . , N 1, where h [n] and hq [n] are the unquantized coefcients respectively and e [n] is the difference between them. This relationship can also be expressed in the frequency domain as Hq () = H () + E () (6.57)
Here H () is the frequency response of the original lter and E () is the error in the frequency response due to coefcient quantization. For IIR lters, the coefcient quantization error may have one more effect: instability. The stability of a lter depends on the location of the roots of the denominator polynomial in the transfer function. Consider a second order section of an IIR lter (since it is the basic building block of higher order lters) with transfer function If the direct form structure is used, assuming rounding, one of the following bounds for the magnitude of the error spectrum can be used: |E ()| = N 2B |E ()| = 2
B
(6.58)
1/2 1/2
(N/3)
(6.59) (6.60)
|E ()| = 2B
1 (N ln N ) 3
where B is the number of bits used for representing the lter coefcients and N is the length of the lter. (6.58) is a worst case absolute bound and is usually overly pessimistic. (6.59) and (6.60) are based on the assumption that the error e [n] is uniformly distributed with zero mean. They generally provide better estimates for the word-length required. Example 6.10. For a length-22 lter with magnitude response shown, at least how many bits are required if the minimum stopband attenuation is 50dB?
160
Figure 6.16: Magnitude Response of an FIR Filter Before and After Coefcient Quantization
161
For IIR lters, the coefcient quantization error may have one more effect: instability. The stability of a lter depends on the location of the roots of the denominator polynomial in the transfer function. Consider a second order section of an IIR lter (since it is the basic building block of higher order lters) with transfer function H (z) = b0 + b1 z 1 + b2 z 2 a0 + a1 z 1 + a2 z 2 (6.61)
The roots of the denominator polynomial, or poles of the transfer function, are located at p1 = p2 1 2 1 = 2 a1 + a1 a2 4a2 1 a2 4a2 1 (6.62) (6.63)
They may either be complex conjugate pairs or are both real. If they are complex conjugate pairs, they can be represented as having a magnitude and an angle: p1 = rej p2 = rej where r = a2 (6.66) (6.67) (6.64) (6.65)
a1 = arccos 2r
For stability, r must be less than 1. This applies to both real and complex poles. So the test for stability for the second coefcient is 0 |a2 | < 1 (6.68)
The arguments to the arc-cosine function in (6.67) must have a magnitude that is less than or equal to 1. So the test for stability for the rst coefcient is |a1 | 2 a2 Both (6.68) and (6.69) must be satised for the IIR lter to remain stable. (6.69)
162
e[n] b0 x[n]
+
z 1 b1 a1 z 1 y[n1] z 1 b2 a2 y[n2] z 1
y[n] + f[n]
x[n1]
x[n2]
Figure 6.17: Second Order IIR Filter Realized in Direct Form I with Round-off Error Included.
6.4.4 Round-off Noise

Round-off noise is introduced when products or sums of products are quantized. Our assumption will be that the rounding off is performed on the result of the (innite precision) sum of product. Consider a second order IIR lter
2 2
y [n] =
k=0
bk x [n k] +
k=1
ak y [n k]
(6.70)
implemented in direct form I. Figure 6.17 shows the the addition of the round-off error modelled by e [n] for this realization. We can view this as a two-input system. The response y [n] is due to input x [n] and the response f [n] is due to input e [n]. Since e [n] is processed only by the right half of the lter
2
f [n] =
k=1
ak f [n k] + e [n]
(6.71)
2 If stationary random white noise e [n] with mean me and variance e is passed through a

e1 [n] e 2 [n]
163
x[n]
+
z 1
+
z 1
y[n] + f[n]
0.25
0.5
Figure 6.18: Cascade Realization of Transfer Function in Example 6.11. system with impulse response h [n], then the output noise f [n] will have a mean
mf = me
n=
h [n] = me H (0)
(6.72)
where H (0) is the DTFT of h [n] at = 0. The variance of the noise is given by
2 f
2 e n=
|h [n]|2
(6.73)
or in the frequency domain

2 f = 2 e 2
|H ()|2 d
(6.74)
Example 6.11. A lter with transfer function H (z) = 1 0.75z 1 1 + 0.075z 2
is to be realized by a cascade of two rst order sections. If the round-off noise variance of the two adders are the same. Find the variance of the round-off noise observed at the 2 output in terms of the round-off noise variance e . Solution:
164
The transfer function can be expressed as a cascade of two rst order sections: H (z) = 1 1 = H1 (z) H2 (z) 1 0.25z 1 1 0.5z 1
Figure 6.18 shows the realization using canonical form. Assuming that the quantization is performed after summation, there are two sources of errors, each introduced by the two adders. The rst error signal e1 [n] passes through the whole system H (z) while the second one e2 [n] only pass through H2 (z). 1 1 1 1 0.5z 1 1 0.25z 2 1 = 1 0.5z 1 1 0.25z 1
n n
H (z) =
So the impulse response is given by h [n] = 2 through an inverse Z transform. Similarly, h2 [n] = 2 The output noise is f [n] = e1 [n] h [n] + e2 [n] h2 [n] The variance of this noise is
2 2 f = e n 2 |h [n]|2 + e n
1 2
u [n]
1 4
u [n]
1 2
u [n]
|h2 [n]|2
2 where e is the variance of the noise signals and is assumed to be the same for both e1 [n] and e2 [n]. Now, 2
|h [n]|
= = 4
2 1 2
1 2
2n
1 4 1 4
n 2
2n
1 8
Finite Wordlength Effects Using the fact that
165
an =
n=0
1 1a
for a < 0, we have 4 1 4 2 + 2 1 (1/8) 1 (0.5) 1 (1/4) = 1.8286
|h [n]|2 =
n
Similarly, |h2 [n]|2 =

n
2 = 2.6667 1 (1/2)2
So
2 2 2 2 f = 1.8286e + 2.6667e = 4.4953e
6.4.5
Limit Cycles
Although we have been treating digital lters as linear systems, the fact is that a digital lter becomes nonlinear when we include quantization, round-off, truncation, and overow. A lter, designed as a linear system, that is stable may oscillate when an overow occurs. This type of oscillation is called a limit cycle. Limit cycles due to round-off, truncation and overow are illustrated by two examples. Example 6.12. A rst order IIR lter has the following different equation: y [n] = x [n] + y [n 1] (6.75)
For = 0.75, the output samples y [n] obtained using initial condition y [0] = 6 and zero input x [n] = 0 for n 0 are listed in Table 6.2. It shows that the output quickly decays to zero. If y [n] is rounded to the nearest integer, then after some time the output remains at 2. The output will now exhibit a DC offset. For = 0.75, the output oscillates briey and decays to zero if computations are done in innite precision. If the output is rounded to the nearest integer, then the output oscillates between 2 and +2. These quantities are listed in Table .
166
n 0 1 2 3 4 5 6 7 8 9 10
y [n], innite precision 6 4.5 3.38 2.53 1.90 1.42 1.07 0.80 0.60 0.45 0.3375
y [n], rounded to nearest integer 6 5 4 3 2 2 2 2 2 2 2
Table 6.2: DC Offset Behaviour of Output Rounding
n 0 1 2 3 4 5 6 7 8 9 10
y [n], innite precision 6 -4.5 3.38 -2.53 1.90 -1.42 1.07 -0.80 0.60 -0.45 0.3375
y [n], rounded to nearest integer 6 -5 4 -3 2 -2 2 -2 2 -2 2
Table 6.3: Oscillating Behaviour of Output Rounding
Finite Wordlength Effects Example 6.13. A stable lter with transfer function H (z) = 1 z 2 z + 0.5
167
is implemented as shown in Figure. Note that we have added an overow non-linearity after the adder in the lter structure. The non-linearity is to model any possible 2s complement overow and the nonlinear function is shown in Figure . If the input is zero, we have x1 [n + 1] = x2 [n] x2 [n + 1] = N L (0.5x1 [n] + x2 [n]) If the initial conditions are x1 [0] = 0.8 x2 [0] = 0.8 Then we have x1 [1] = 0.8 x2 [n] = N L (0.5x1 [0] + x2 [0]) = N L (1.2) = 0.8 In fact, for n 1, x1 [n] = (1)n 0.8 x2 [n] = (1)n+1 0.8 Thus the output oscillates between 0.8 and 0.8. Note that the limit cycle will only start if there is a previous overow. If no overow occurs, the system remains linear. (6.80) (6.81) (6.78) (6.76) (6.77)
(6.79)
6.4.6
Arithmetic Overow
Overow occurs when the two large number of the same sign are added and the result exceeds the word-length. If we use 2s complement representation, as long as the nal
168
x[n]
NL
z 1
x 1 [n]
z 1
x 2 [n]
y[n]
0.5
+
Figure 6.19: Filter Structure with Overow Nonlinearity. result is within the word-length, overow of partial results is unimportant. If the nal result does cause overow, it may lead to serious errors in the system. Overow can be avoided by detecting and correcting the error when it occurs. However, this is a rather expensive approach. A better way is to try to avoid it by scaling the data or the lter coefcients. Overow can potentially cause an IIR system to become unstable. It may also lead to limit cycles similar to round-off errors. Consider a lter with transfer function H (z) = 1 z 3 z + 0.5 (6.82)
It can easily be veried that it is a stable system. A realization of this lter is shown in Figure 6.19. Here we have added a nonlinear function N L () to the realization to model any potential 2s complement overow that may occur as a result of the addition. This nonlinear function is shown in Figure 6.20. Consider the case where the input is zero with the initial conditions x1 [0] = 0.8 x2 [0] = 0.8 Then from Figure 6.19 we have, for n 0, x1 [n + 1] = x2 [n] x2 [n + 1] = N L (0.5x1 [n] + x2 [n]) (6.85) (6.86) (6.83) (6.84)
Summary
NL(x)
169
Figure 6.20: 2complement Overow Modelled as a Nonlinear Function. So x1 [1] = 0.8 x2 [1] = 0.8 and in general, for n 1, x1 [n] = (1)n 0.8 x2 [n] = (1)n+1 0.8 Thus the lter is stuck in a limit cycle. (6.89) (6.90) (6.87) (6.88)
6.5
Summary
The basic IIR lter equation (which includes the FIR equation as a special case) looks simple. However, there are a number of ways by which it can be realized in both hardware and software. Each realization has its chacteristics. The basic forms are the direct forms I and II and their transposed counterparts. A higher order lter is usually broken up into rst and second order sub-systems organized in a cascade or parallel form. In digital signal processors, apart from optimizing the MAC operation, the use of circular buffers
170
also has a signicant impact on ltering. We have also looked at a number of practical implementation issues including nite wordlength effects on the lter frequency response and the operational issues like overow. Whether the computations are implemented in VLSI/ASIC/FPGA or in software running on digital signal processors or general purpose processors, the issues are still relevant. In the following chapters, we shall be concerned with the design of FIR and IIR lters.
Chapter 7 FIR Filter Design

Finite Impulse Response (FIR) Filters are digital lters with impulse responses that vanish after a nite number of samples. The output of the difference equation of an FIR lter depends only on a nite number of past and present inputs. The lter design problem is to determine the lter coefcients that produce a lter transfer function which meets the prescribed frequency response specications. For FIR lters, this is equivalent to designing the nite impulse response sequence of the lter. The main advantage of FIR lters is that it can have exact linear phase. We shall assume that only linear phase lters are required.
7.1
Design Process
The general digital lter design process can be broken up into four main steps: 1. Approximation 2. Synthesis and realization 3. Performance analysis 4. Implementation The process is illustrated in Figure 7.1.
172
FIR Filter Design
Filter Specifications
Approximation Filter Transfer Function Synthesis and Realization Filter Structures
Choose alternative structure if necessary
Performance Analysis Stability Quantization error Implementation
Filter implemented as hardware and/or software
Figure 7.1: Digital Filter Design Process
Design Process
173
Specications The design process normally starts with the specications and requirements of the lter which are intimately related to the application at hand. These specications may include frequency domain characteristics such as magnitude and phase responses. There may also be some time domain requirements such as maximum delay. Most specications dene the upper and lower limits to each of these characteristics either in absolute terms or in decibels. They are illustrated in Figure 7.2. Alternatively, a desired or ideal response may be given with the maximum amount of deviations from the ideal specied. Passband ripple in dB is given by Rp = 20 log10 and the stopband attenuation is given by As = 20 log10 s 1 + p (7.2) 1 p 1 + p (7.1)
Approximation Given the lter specications, the rst step of the design process is to nd a lter transfer function that will satisfy these specications. This process is called approximation. It is so called because what we are doing is in fact nding a transfer function that approximates the ideal response that is specied. The methods for solving the approximation problem for digital lters can be classied as direct or indirect. With direct methods, the problem is solved in the discrete-time (and hence discrete-frequency) domain. For indirect methods, a continuous-time transfer function is rst obtained using well-established methods in analog lter design. This transfer function is then transformed into a discrete-time transfer function. Indirect methods are more commonly used for IIR lters whereas FIR lter design methods are mostly direct ones. These solution methods can also be classied as closed-form or iterative. Closed form methods make use of closed-form formulas and are usually completed in a denite number of steps. Iterative methods make use of optimization techniques that starts with an initial solution and the solution is rened progressively until some pre-determined performance criteria are satised. The number of iterations are unknown and depends on the initial solution and the effectiveness of the optimization techniques employed. The following methods will be discussed in more detail in Sections 7.2 and 7.5:
174
FIR Filter Design
|H( )| passband ripple allowed
1+ p 1 p
transition band
s 0 0
p
minimum stopband attenuation
(a) Specications in Absolute Values
|H( )| (dB) passband ripple allowed
0 Rp
transition band
As 0
p
minimum stopband attenuation
(b) Specciations in dB
Figure 7.2: Filter Specications in Frequency Domain.
Design Process 1. Window method. 2. Frequency sampling method. 3. Parks-McClelland (Remez Exchange) algorithm.
175
Realization Once the transfer function has been determined, it has to be realized in some form. This procedure is analogous to the lter realization procedure for analog lters where a suitable circuit topology and circuit element values are chosen to realize a certain lter transfer function. A number of realization methods has been discussed in Chapter 6. The best realization of a given transfer function depends very much on the application. General considerations include the number of adders and multipliers required, and the sensitivity of the network to nite precision arithmetic effects.
Performance Analysis Even though the lter coefcients are determined to a high degree of precision in the approximation step, digital hardware has a nite precision. The accuracy of the output will depend on the type of arithmetic used: xed-point or oatingpoint. This is particularly so for xed point arithmetics. The designer must ensure that the error introduced by nite precision will not cause violations of the lter specications. Furthermore, arithmetic overow and underow effects must be examined. It cannot be over-emphasized how important this design step is, especially for IIR lters. While FIR lters are guaranteed to be stable, IIR lters can exhibit instability due to quantization errors introduced in the computational process.
Implementation Digital lters can be implemented either in software or hardware or a combination of both. Software implementations require a decision to be made on the type of computer or microprocessor the software will eventually run on. DSP chips which are designed specically for DSP type of operations are very effective. In very demanding applications, the lter may need to be hard-wired or implemented as an application specic integrated circuit (ASIC) in order to obtain the speed required. It may also be necessary that some of the other functions such as analog-to-digital conversion and digital-to-analog conversion be integrated on the same device. However, development time will generally be longer and the cost is much higher.
176
FIR Filter Design
7.2
Window Method
The window method is suitable for designing lters with standard low-pass, bandpass and high-pass responses. The method outline is as follows: 1. Choose an ideal lter response. 2. Calculate its (innite) impulse response using inverse discrete-time Fourier Transform. 3. Multiply the impulse response with a nite-length window sequence. 4. If a causal lter is required, time shift this nite-length impulse response so that all the non-zero values are in the non-negative part of the time axis. Consider an ideal lowpass lter with cutoff frequency c . The frequency response of this ideal lter is ej , || c 0, c < || <
Hd () =
(7.3)
Notice that we have introduced an arbitrary phase shift in the frequency response. Using inverse DTFT, we obtain the impulse reponse of this ideal lter: hd [n] = sin [c (n )] (n ) (7.4)
This impulse response has innite length. In order to obtain a lter with length N , hd [n] is windowed (or multiplied) by a window of length N . The impulse response h [n] of the lter is zero for n < 0 and n N . Some of the window functions that can be used are essentially the same as those used for spectral analysis. The rectangular and Hamming windows have been discussed in Section 4.3. Other ones include the Hanning window, Blackman window and Kaiser window. The phase shift in the frequency response corresponds to a delay in the impulse response. This is useful if we want to design lters that are causal.
Window Method
177
7.2.1 Rectangular Window

The second step in the approximation process using the windowing method is to window the innitely long ideal impulse reponse with a nite length window. We shall start with the rectangular window due to its simplicity. The following example illustrates the process. Example 7.1. Design a 10-th order linear phase FIR lter with cutoff frequency of /4 radians using a rectangular window. Solution: In this case, N = 10 and so = N =5 2
The windowed impulse response is given by (7.4) and hd [n] = sin c (n ) 0nN (n ) sin (n 5) 4 = (n 5) 2 2 1 2 1 2 1 2 2 = , 0, , , , , , , , 0, 10 6 2 2 4 2 2 6 10
How good is the design using rectangular windows? In other words, how closely does the resulting FIR lter frequency response H () approximate the original ideal response Hd ()? To answer this question, we performed the DTFT of hd [n] obtained in Example 7.1 and plotted the magnitude and phase responses in Figure 7.3. Note that we indeed have exact linear phase for this lter since the lter coefcients are symmetric. Ripples exist in the stopand and the minimum stopband attenuation is about 20dB. Figure 7.4 shows the FIR lter magnitude responses for N = 51 and N = 101. Notice that there are now some ripples in the passband also. One would expect that as N increases, the approximation will become better. This is indeed the case except for the region close to the transition between passband and stopband. In this region the Gibbs phenomenon discussed in Section 3.3.1 manifests itself. For this reason, the approximation at the band edge will always be poor for the rectangular window design regardless of how large N is.
178
10th Order FIR Filter Response
FIR Filter Design
0 10 Magnitude (dB) 20 30 40 50 60 70 0 0.1 0.2
0.3 0.4 0.5 0.6 0.7 Normalized Frequency ( rad/sample)
0.8
0.9
0 50 Phase (degrees) 100 150 200 250 300 350 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Normalized Frequency ( rad/sample) 0.8 0.9 1
Figure 7.3: FIR Filter Response for N = 11. We can now summarize the characteristics of using the rectangular window: The ripple size decreases with increasing lter order N . Approximation well within the passband and stopband becomes better as the lter order increases. The transition width decreases with increasing lter order. For any order N , the lter response is always equal to 0.5 at the cutoff frequency. Ripple size near the passband to stopband transition remains roughly the same as N increases. The maximum ripple size is about 9%. This is known as the Gibbs phenomenon. Using the Hamming window will provided better passband responses, especially at the band edges. Other window functions that can be used include
Window Method
179
20 0 Magnitude (dB) 20 40 60 80 100 0 0.1 0.2
0.8
0.9
0 200 Phase (degrees) 400 600 800 1000 1200 1400 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Normalized Frequency ( rad/sample) 0.8 0.9 1
(a) 50th Order FIR Filter Responses
20 0 Magnitude (dB) 20 40 60 80 100 0 0.1 0.2
0.8
0.9
0 500 Phase (degrees) 1000 1500 2000 2500
0.1
0.2
0.8
0.9
(b) 100th Order Filter Responses
Figure 7.4: FIR Filter Responses for N = 51 and N = 101.
180 Hanning window w [n] = Blackman window w n= 0.42 0.5 cos 0,

2n N 1
FIR Filter Design
0.5 1 cos 0,
2n N 1
, 0nN 1 otherwise
(7.5)
+ 0.08 cos
4n N 1
, 0nN 1 otherwise
(7.6)
These two windows are variants of the Hamming window; they both possess the raised cosine shape. The transition band width and the stopband attenuation of the resulting lter is largely dependent on the window that is used. In order to design lters with different parameters, a unied window function would be useful. The Kaiser window is one such window function.
7.2.2
Kaiser Window
The Kaiser window is the most commonly used window when it comes to designing FIR lters using the window method. It is actually a class of window function that is parametrized by a variable . It is dened as I0 x [n] = 1 1 I0 ()
2n 2 N 1
0nN 1
(7.7)
where I0 () is the modied zero-order Bessel function dened by the series expansion
I0 (x) = 1 +
k=1
x k 2
k!
(7.8)
This window is optimal in the sense that it provides the sharpest transition band width with the largest main lobe. can be chosen to give various transition widths and stopband attenuation. For example, if = 6.658, then the transition width is 7.8/N and the minimum stopband attenuation is about 60 dB. When we set = 4.538, the transition width becomes 5.8/N with a
Highpass Filters
181
minimum stopband attenuation of around 50 dB. One accepted way of determining is to use one of the following formulas, depending on the stopband attenuation As : As 21 0, 0.4 = (7.9) 0.5842 (As 21) + 0.07886 (As 21) , 21 < As 50 0.1102 (As 8.7) As > 50 The lter order can be determined by N 1 0 D (s p ) (7.10)
where p and s are the passband and stopband edge frequencies and 0 is the sampling frequency. The parameter D is given by D= 0.9222,
(As 7.95) , 14.36
As 21 As > 21
(7.11)
7.3 Highpass Filters

High-pass lter design using Kaiser windows is very similar to low-pass lter design. The only change in the steps is simply dene f = fpass fstop since fpass > fstop for high-pass lters. The ideal high-pass impulse response is obtained from the inverse DTFT of the ideal high-pass frequency response. It is given by hd [k] = [k] sin (c k) k (7.13) (7.12)
The windowed lter impulse response is therefore h [n] = w [n] [n M ] sin ((n M ) c ) (n M ) sin ((n M ) c ) = [n M ] w n (n M )
(7.14)
182
FIR Filter Design

lowpass output
HLP (z)
highpass output
Figure 7.5: Highpass Filter Implemented by a Lowpass Filter and Delay. The second term on the right hand side of this equation is the impulse response of the low-pass lter with the same cutoff frequency. Note that with the same value of c , the low-pass and high-pass lters are complementary. That is, hLP [n] + hHP [n] = [n M ] In terms of Z-transform, (7.15) becomes HLP (z) + HHP (z) = z M (7.16) n = 0, 1, . . . , N 1 (7.15)
In other words, we can obtain a high-pass lter by designing a low-pass lter with the same cutoff frequency and using (7.15) above. Alternatively, the high-pass lter can be implemented using a low-pass lter and a delay as shown in Figure .
7.4
Bandpass Filters
Bandpass lters can also be designed via low-pass lters. The ideal bandpass lter frequency response is given by HBP () = 1, 2p || 1p 0, otherwise (7.17)
This response can be expressed as the difference between two low-pass responses HBP () = HLP 1 () HLP 2 () (7.18)
Frequency Sampling Method
183
where HLP 1 () is a low-pass response with cutoff frequency 1p and HLP 2 () is one with cutoff frequency 2p . These low-pass lters have impulse responses hLP 1 [n] = hLP 2 [n] =
sin 1p n , n 1p , sin 2p n , n 2p ,
n=0 n=0 n=0 n=0
(7.19) (7.20)
The bandpass lter impulse response is therefore given by hBP [n] =

sin 2p nsin 1p n , n 2p 1p ,
n=0 n=0
(7.21)
7.5 Frequency Sampling Method

Another relatively simple FIR lter design method is the frequency sampling method. While the windowing method is only suitable for designing lters with simple standard frequency responses like low-pass, high-pass and band-pass, the frequency sampling method can handle arbitrary frequency responses. The basic idea is that the nite impulse response of an FIR lter can be obtained from the DFT of its frequency response. So if we sample the frequency response at N regularly spaced frequency points (between to ), then we can use these frequency samples to calculate the impulse response. If samples of the desired frequency response is given as Hd [2k/N ], then the lter coefcients h [n] can be obtained by solving the following N simultaneous equations:
N 1
h [n] ej2nk/N = Hd
n=0
2k N
(7.22)
The coefcients can also be computed via the IDFT which is straight forward. The main problem of the frequency sampling method lies in the determination of Hd [2k/N ] since they are complex-valued in general. Let Hd 2k = A [k] ej[k] N (7.23)
where A [k] and [k] are samples of the amplitude and phase response respectively. If we want the lter to have linear phase, then we need to have symmetrical or anti-symmetric h [n].
184 For symmetrical h [n] and the order of the lter N odd, we require A [k] = A [N k] [k] = k (N 1) N 1k N 1 2
FIR Filter Design
0k<N
In this case, the impulse response is given by (N 1)/2 2 k (1 + 2n) h [n] = A [0] + 2 (1)k A [k] cos N 1 N k=1 for 0 n N 1. If N is even, then A N 2 = 0 N 2 k(N 1) N , 0k N 2 k(N 1) , N < k < N N 2 1k
(7.24)
(7.25) (7.26) (7.27)
A [k] = A [N k] [k] = The impulse response is h [n] = for 0 n N 1. 2 A [0] + 2 N 1

N 2
(1)k A [k] cos k (1 + 2n) N (7.28)
k=1
While the frequency sampling method is more exible compared to the windowing method, care should be taken with the transition region to make sure that there are enough frequency samples to constrain the transition band behaviour.
7.6 Summary
The FIR digital lter design process is outlined in this chapter. Two design (or approximation) methods are described: windowing method and frequency sampling method. The
Summary
185
windowing method is suitable for designing standard low-pass, high-pass and bandpass lters. The frequency sampling method is more exible in that any arbitrary frequency response could be specied. However, care must be taken to ensure that the transition band behaviour conforms to what is expected. The interested reader is encourage to consult the vast number of books available on digital lter design on other design methods.
186
FIR Filter Design
Chapter 8 IIR Filter Design

While FIR lters have impulse responses that vanish after a nite number of samples. Innite Impulse Response (IIR) Filters are ones with the output sample depending on past input and past output samples.
N M
y[n] =
k=1
a[k]y[n k] +
k=0
b[k]x[n k]
(8.1)
The rst term on the right-hand-side of the difference equation (8.1) is the feedback of past outputs to the current output. It is this recursive part of the equation that causes the impulse response to be inntely long. For this reason, IIR lters are also called recursive lters. Since the IIR lter is a linear system, the output and input are related by the convolution sum y [n] =
k=0
h [k] x [n k]
(8.2)
assuming a causal lter. What makes it different from FIR lters is that the upper limit of the summation is innity because the impulse response h [n] has innite duration. Obviously computation using (8.2) is impractical. The recursive relationship dened in (8.1) is much more efcient and practical. It requires N + M + 1 multiplications and N + M additions if implemented directly. Since the impulse response is innitely long, we cannot adapt the techniques for FIR lter design to IIR lters. The most common techniques involve rst designing an analog lter with equivalent specications.
188
Digital Filter Specifications Analog Filter Specifications
IIR Filter Design
Conversion
Analog Filter Design Methods Digital IIR Filter Analog Filter
Transformation
Figure 8.1: IIR Filter Design Approach
8.1
Design Approach
One of the most efcient ways of designing IIR digital lters is through a corresponding analog lter. The process is illustrated in Figure 8.1. First, the digital lter specications are translated to analog lter specications. Then an analog lter is designed according to those specications. The transfer function of the resulting analog lter is then transformed into an equivalent digital lter transfer function. We have a wealth of knowledge about analog lter approximations which we can make use of. The only step that requires some thought is the transformation from an analog lter to a digital one. In this book, only two methods are discussed. They are the impulse invariance method and the bilinear transformation method.
8.2
Some Classic Analog Filters
We shall briey review two types of classic analog lters: Butterworth and Chebyshev. The interested reader is encouraged to consult textbooks on other types of analog lters.
8.2.1
Butterworth Filters
Butterworth lters have magnitude squared frequency response given by
189
|H ()|2 =
1 1 + (/0 )2n
(8.3)
where n (an integer) is the order of the lter and 0 is some constant. Designing a Butterworth lter is equivalent to computing the values of n and 0 for a given specication. Consider a low-pass Butterworth lter with passband edge at p and stopband edge at s . Suppose the passband ripple Kp and stopband attenuation Ks are given in dB, then for = s 1 |H (s )|2 = (8.4) 1 + (s /0 )2n and so Ks = 10 log10 |H (s )| = 10 log10 1 + which gives Ks = log10 1 + 10 Similarly for the = p , Kp = 10 log10 |H (p )|2 = 10 log10 1 + and therefore Kp = log10 1 + 10 p 0
2n 2
s 0
2n
(8.5)
s 0
2n
(8.6) (8.7)
s 0
2n
= 100.1Ks 1
(8.8)
2n
p 0
(8.9)
p 0
2n
(8.10) (8.11)
= 100.1Kp 1
190 Combining (8.7) and (8.11), we obtain s p s 2n log p 100.1Ks 1 100.1Kp 1 100.1Ks 1 = log 100.1Kp 1 log M n = log = s p 100.1Ks 1 100.1Kp 1
2n
IIR Filter Design
(8.12) (8.13) (8.14)
where = M = (8.15) (8.16)
The order of the lter is the integer closest to but higher than n obtained from (8.14). Example 8.1. Determine the order of a Butterworth lter that meets the following specications: p s Ks Kp Solution: = 1.5 M = 103 1 = 793.534 100.1 1 log10 793.534 n = 16.466 log10 1.5 = = = = 1 rad/s 1.5 rad/s 30dB 1dB
So the minimum order of the Butterworth lter that satises these specications is 17. Figure 8.2 shows the magnitude response of a 17-th order Butterworth lter. Once the order of the lter is determined, how do we obtain the transfer function H (s) of the lter? Notice that given n, only |H ()|2 is known.

50 Magnitude Response of a Normalized 17th Order Butterworth Filter
191
Magnitude (dB)
50
100
150
200
0.5
1.5 2 2.5 Angular Frequency (normalized)
3.5
Figure 8.2: Magnitude Response of a 17-th Order Butterworth Filter However, we know that the poles of |H ()|2 are located in the same place as H (s) H (s). So H (s) H (s) = 1+ = and s = [ (1)n ]
1/2n
1
0 2n
0 2
(8.17)
=s2
1 1 + (1)n s2n
(8.18)
(8.19)
They are evenly spread out along the unit circle. The angles of the poles will be the same on the s-plane for H (s) H (s) while the distance from the origin is scaled by 0 . For the lter to be a stable system, all its poles have to be on the left hand side of the s-plane. Therefore H (s) will only consists of those poles on the left hand side of the imaginary axis. Thus it can be computed.
192
IIR Filter Design
Example 8.2. Design a Butterworth low-pass lter with the following specications: p s Ks Kp Solution: = 2 M = 102 1 = 19.554 100.1 1 log10 19.5354 = 4.289 n log10 2 = = = = 1 rad/s 2 rad/s 20dB 1dB
So we need a 5-th order lter. At = p p 0 according to (8.11) and 10 (log p log 0 ) = log 0.2589 log 0 = log 1 0.1 log 0.2589 = 0.05868 0 = 1.145 The poles of the normalized Butterworth lter are located as shown in Figure 8.3.
10
= 100.1 1
(8.20)
(8.21)
8.2.2
Chebyshev Filters
The Butterworth lter, while possessing some desirable properties, does not provide a sufciently good approximation near the passband edge. The roll-off from passband to stopband is also relatively gradual. Hence for lters that have a narrow transition band, a very high order Butterworth lter is required. If the application can tolerate some ripples in the passband, the Chebyshev approximation will overcome some of the problems
193
1 0.8 0.6 0.4 0.2 0 0.2 0.4 0.6 0.8 1
poles
0.8
0.6
0.4
0.2
0.2
0.4
0.6
0.8
Figure 8.3: Pole Locations for H(s)H(-s) of a Normalized 5th Order Butterworth Filter.
194
IIR Filter Design
associated with Butterworth lters stated above. The magnitude squared response of Chebyshev lters has the form |H ()|2 = 1 1+
2C 2 N
()
(8.22)
where CN () is the N -th order Chebyshev polynomial and is a parameter associated with the size of the ripples. The Chebyshev polynomials can be obtained recursively: C1 () = C2 () = 2 2 1 Cn+1 () = 2Cn () Cn1 () These polynomials have the following properties
2 1. 0 Cn () 1 for 0 1. 2 2. Cn () 1 for 1.
(8.23) (8.24) (8.25)
That is, they oscillate between 1 and +1 for 1 < < 1 and increase monotonically outside this interval. So the low-pass Chebyshev lter response has a normalized passband from 0 to 1 rad/s. Figure 8.4 shows the magnitude response of a 7-th order Chebyshev lter. Designing a Chebyshev lter involves determining both the order N and the parameter . Generally the order inuences the transition width and the number of oscillations within the passband and inuences the size of the ripples in the passband. Within the normalized passband, 1 1+
2
|H ()|2 1
(8.26)
Thus the specied maximum allowable passband ripple Kp of the lter is related to by Kp = 10 log10 1 + or = 100.1Kp 1 (8.28)
2
(8.27)

10 Magnitude Response of a 7th Order Chebyshev Filter
195
10
20 Magnitude (dB)
30
40
50
60
70
80
0.5
1 1.5 Angular Frequency (normalized)
2.5
Figure 8.4: Magnitude Response of a 7-th Order Chebyshev Filter with 2dB Passband Ripple. Without going into the details of the derivation, the order of the lter given a specied minimum stopband attenuation of Ks dB with passband edge at p and stopband edge at s is given by cosh1 M N= cosh1 where = M = similar to that for Butterworth lters. s p 100.1Ks 1 100.1Kp 1 (8.30) (8.31) (8.29)
196
IIR Filter Design
8.3
Bilinear Transformation
Once the analog lter has been designed (i.e. H (s) has been determined), we need to map it to the digital domain to obtain H (z). Bilinear transformation is one such method for converting the analog lter transfer function H(s) into a corresponding digital one H(z). The transformation is performed by a change of variables s= 1 + z 1 1 z 1 (8.32)
It is called bilinear because both the numerator and denominator of this transformation equation are linear. This transformation is reversible in that H (s) can be obtained from H(z) by the substitution 1+s z= (8.33) 1s Through the bilinear transformation, the imaginary axis of the s-plane is mapped to the unit circle of the z-plane. This can be shown by letting s = + j and so through (8.33) we get 1 + + j (8.34) z= 1 j If = 0, i.e. s = j (the imaginary axis of the s-plane), then z= 1 + j 1 j (8.35)
Obviously the magnitude of the numerator and denominator are the same and so |z| = 1 which is the unit circle on the z-plane. If we consider < 0 in (8.34), then (1 ) > (1 + ) and |z| < 1. So the bilinear transformation also maps the left hand side of the s-plane to the inside of the unit circle on the z-plane. (8.36)
8.3.1
Frequency Warping
In order to understand the effects of this transformation, let and be the frequency variables in the analog (s) and digital (z) domains respectively. The frequency response
Bilinear Transformation of the digital lter can be obtained from H(z) by the substitution z = ejt and that of the analog lter by substituting s = j into H (s). Using the above substitutions in (8.32), we have 1 ej j = 1 + ej ej/2 ej/2 ej/2 = j/2 j/2 e (e + ej/2 ) 1 ej/2 ej/2 2j = j 1 j/2 (e + ej/2 ) 2 which can be simplied to = sin (/2) cos (/2) = tan 2
197
(8.37)
(8.38)
(8.39)
(8.40)
This relationship is plotted in Figure 8.5. It shows three bands of frequencies in the analog domain that are of equal width. After applying the bilinear transformation, the corresponding bands of frequencies in the digital domain are not of the same width. This nonlinear relationship between the analog and digital frequency variables lead to a distortion or warping of the digital frequency response. This means that the passband of an analog lter will not be the same as that of the digital lter obtained through bilinear transformation. Given the passband specications of the digital lter, the passband frequencies of the correct analog lter is obtained by pre-warping. Suppose fc is a critical frequency (e.g. cutoff frequency, passband edge, stopband edge, etc.), then the corresponding pre-warped analog frequency is obtained by c = tan fc fs (8.41)
where fs is the sample frequency of the digital lter. Let us illustrate the design process through an example.
198
Analog Frequency
IIR Filter Design
Equal Width Digital Frequency
Unequal Width
Figure 8.5: Frequency Warping in Bilinear Transformation Example 8.3. Design a low-pass IIR digital lter with the following specications: 0.89125 |H| 1, |H ()| 0.17783, 0 || 0.2 0.3 ||
Characteristics of the bilinear transformation are summarized below: 1. Provided that the analog lter is stable, the resulting digital lter is guaranteed to be stable. 2. The order of the digital lter is the same as the prototype analog lter. 3. Optimal approximations to the analog lter transform into optimal digital lters.
Bilinear Transformation
199
4. The cascade of sections designed by the bilinear transformation is the same as that obtained by transforming the total system.
8.3.2
Frequency Transformation
So far we have only considered low-pass lter designs. Designing high-pass and bandpass lters requires a transformation to be performed on the low-pass lter. Three approaches can be used for IIR digital lters. 1. In the rst approach, the analog prototype lter is transformed to be appropriate high-pass or band-pass analog lter. Then this analog lter is turned into an IIR lter by bilinear transformation. The low-pass to high-pass transformation is given by 1 s (8.42) s and the low-pass ot band-pass transformation is s s2 + u l s (u l ) (8.43)
where the low-pass lter has a cutoff frequency of 1 rad/s and u and l are the upper and lower passband edges of the bandpass lter respectively. 2. Alternatively, the analog low-pass lter is rst converted to a digital low-pass lter. The digital high-pass or band-pass lter is then obtained by a spectral transformation of the digital low-pass lter. This method is less commonly used. 3. The last approach is basically a combination of the two previous methods. New bilinear transform equations can be obtained that directly maps an analog low-pass lter to a digital high-pass, bandpass or band-reject lters. The bilinear transform function for low-pass to high-pass is s= 1 + z 1 1 z 1 (8.44)
and the low-pass to bandpass transform is s= where c is a certain constant. 1 2cz 1 + z 2 1 z 2 (8.45)
200
IIR Filter Design
Example 8.4. Design a high-pass IIR lter operating at a rate of 10 kHz. The stopband and passband edges are 2 kHz and 4 kHz respectively. The maximum passband attenuation is 0.5 dB and the minimum stopband attenuation is 10 dB. Use Butterworth lter as the analog prototype lter and the low-pass to high-pass bilinear transformation 1 + z 1 s= 1 z 1 Solution: First establish the frequency warping function between the analog and digital frequencies. Substituting s = j and z = ej into the above bilinear function, we have j = 1 + ej 1 ej 1 1 ej/2 ej/2 + ej/2 2 = 1 j 2j ej/2 (ej/2 ej/2 )
1 cos (/2) j sin (/2) = cot 2 = Now, fs = 10kHz fpass = 4kHz p = fstop 2 (4) = 0.8 10 2 (2) = 2kHz s = = 0.4 10
(8.46)
Using (8.46) we can obtain the corresponding analog band edge frequencies p = cot s = cot 0.8 2 0.4 2 = 0.3249 = 1.3764
Since the magnitude spectrum of the lter is symmetrical about f = 0, we do not need to worry about the negative signs. Hence the passband edge frequency of the analog (low-pass) lter is 0.3249 radians and the stopband edge frequency is 1.3764 radians.
Bilinear Transformation The Butterworth lter response is given by |H ()|2 = For = p , and for = s , 1 1 + (/0 )2n
201
(8.47)
|H (p )|2 = 100.5/10 = 100.05 |H (s )|2 = 1010/10 = 0.1
according to the specications for passband and stopband attenuations. Substituting these values into (8.47) gives p = 100.05 1 0 2n (log p log 0 ) = log 100.05 1 2n (0.4883 log 0 ) = 0.9136 and s = 10 1 = 9 0 2n (log s log 0 ) = 0.9542 0.9542 log 0 = 0.1387 2n Substituting (8.49) into (8.48), we get 2n 0.4883 0.1387 + 0.9542 = 0.9136 2n (0.1254n + 0.9542) = 0.9136 n = 1.48
2n 2n
(8.48)
(8.49)
So the order the Butterworth lter required is 2. Now we can compute the value of 0 . log 0 = 0.1387 0 = 0.7946 0.9542 = 0.09985 2 (2)
202
IIR Filter Design
The poles of the Butterworth lter are those on the left-hand-side of the s-plane, i.e. at cos and so H (s) =
s 0 2
j sin 4 4
= 0.7071 j0.7071 1 + 2
(8.50)
s 0
+1
To obtain the digital IIR lter transfer function, we substitute s= into (8.50): H (z) = 1.5838 = 1
1+z 1 1z 1 2 1+z 1 + 1z 1 1 2
1 + z 1 1 z 1
+ 1.778
(1 z ) 1.5838 (1 + + 1.778 (1 + z 1 ) (1 z 1 ) + (1 z 1 )2 1 2z 1 + z 2 = 1.5838 + 3.1676z 1 + 1.5838z 2 + 1.778 1.778z 2 + 1 2z 1 + z 2 1 2z 1 + z 2 = 4.3618 + 1.1676z 1 + 0.8058z 2 (0.2293) (1 2z 1 + z 2 ) = 1 + 0.2677z 1 + 0.1847z 2 z 1 )2
8.4
Impulse Invariant Method
Let Ha () be the transfer function of the analog lter that has been designed. The impulse response ha (t) of this lter can be obtained through Fourier transformation. The idea behind the impulse invariance method is that the impulse response of the digital lter h [n] is a sampled version of ha (t). Thus h [n] = ha (nT ) where T is the sampling period. (8.51)
203
As we know, the magnitude spectrum of a sampled signal is a periodic continuation of the original spectrum of the analog signal. So the magnitude response H () of the digital lter is periodic with a period of fs = 1/T . It can be expressed mathematically as 1 H () = T
Ha
k=
j ( 2k) T
(8.52)
Even though the impulse invariant method preserves the shape of the impulse response, the frequency response may be signicantly different from what we expected. This is because the stopband magnitude response does not drop to zero at = . This means that H() will be an aliased version of Ha () because of the periodic nature of the response. This is illustrated in Figure 8.6. Due to the aliasing problem, the impulse invariant technique is only applicable to low-pass lter designs.
|()|
Analog Filter Response
Replica of Analog Filter Response
Combined Response
Figure 8.6: A Sampled Magnitude Response
The stopband attenuation of the aliased version of Ha () may not be sufcient to satisfy the original lter specications. This will render the design useless. So we have to make sure that the aliased portion of H() is small enough. The passband is also affected but the effect is usually much smaller than that in the stopband. A sufciently high sampling
204
IIR Filter Design
rate can be used to overcome this problem; the higher the sampling rate, the smaller the error introduced by aliasing. For this reason, Butterworth and Chebyshev lters are more suitable if impulse invariant method is used. This is because both these lters do not have ripples in the stopband and the response in this region is monotonically decreasing. The analog lter transfer function can be expressed in terms of partial fractions
N
Ha () =
i=1
Ai j si
(8.53)
where si is generally complex-valued. In fact, it is usually expressed in terms of the Laplace variable s: N Ai Ha (s) = (8.54) s si i=1 Inverse Laplace transform of (8.54) gives the impulse response ha (t) =
N i=1
0,
Ak esk t , t 0 t<0
(8.55)
The impulse response of the digital lter is obtained by sampling. Therefore, h [n] = ha (t)
N
(8.56) (8.57) (8.58)
=
i=1 N
Ak esk nT u [n] Ak esk T

i=1 n
= Applying the Z transform pair
u [n]
an u [n] we obtain H (z) =

i=1 N
1 1 az 1 Ai 1 esk T z 1
(8.59)
(8.60)
205
Since the transfer function of the digital lter can be obtained directly from Ha (s), the coefcients of the IIR lter can be obtained directly from Ha (s) without having to compute the impulse response rst. In summary, the characteristics of impulse invariance designed lters are: 1. The IIR lter is stable if the analog lter is stable. 2. The frequency response of the IIR lter is an aliased version of the analog lter. So the optimal properties of the analog lter are not preserved. 3. The cascade of two impulse invariance designed lters is not impulse invariant with the cascade of the two analog lter prototypes. 4. The step response (the response of the lter to a unit step input) is not a sampled version of the step response of the analog lter. Therefore, if the step response of the nal lter is important, then the impulse invariant method may not be a suitable choice. Example 8.5. A second order Butterworth lter has the following normalized transfer function 1 H (s) = 2+ s 2s + 1 Design an IIR lter based on this analog lter using the impulse invariant method. Solution: First expand the transfer function in terms of partial fractions: H (s) = s1 = s2 = A1 = A2 = A1 A2 + s s1 s s2 1 (1 j) 2 1 (1 + j) 2 j 2 j 2
206 Then the IIR lter transfer function is given by H (z) = A1 A2 + s1 T z 1 1e 1 es2 T z 1
IIR Filter Design
Since this is a normalized transfer function, we can use the normalized sampling period T = 1. So the transfer function becomes H (z) = 0.2265z 1 1 0.7497z 1 + 0.2431z 2
8.5 Summary
Digital IIR lters designs are usually based on an analog lter with equivalent specications. We have briey reviewed two common types of analog lters: Butterworth and Chebyshev. Two methods for converting the analog lter response to a digital one has been discussed. The impulse invariant method seeks to preserve the shape of the impulse response of the analog lter. However, it is only suitable for designing low-pass lters because of aliasing introduced when sampling. The second method, bilinear transformation, can be applied to a larger class of lters. However, the frequency warping effect has to be compensated for.
Bibliography
[1] S. J. Orfanidis, Introduction to Signal Processing, ser. Prentice-Hall Signal Processing Series. Prentice-Hall, 1996. [2] A. Oppenheim and R. Schafer, Discrete-Time Signal Processing, second edition ed., ser. Prentice-Hall Signal Processing Series. Prentice-Hall, 1999. [3] S. K. Mitra, Digital Signal Processing: A Computer-Based Approach, second edition ed. McGraw-Hill, 2001. [4] J. G. Proakis and D. G. Manolakis, Digital Signal Processing: Principles, Algorithms and Applications, third edition ed. Prentice-Hall, 1996. [5] C. Marven and G. Ewers, A Simple Approach to Digital Signal Processing. WileyInterscience, 1996. [6] R. G. Lyons, Understanding Digital Signal Processing. Addison-Wesley, 1997.
[7] J. Y. Stein, Digital Signal Processing: A Computer Science Perspective, ser. Wiley Series in Telecommunications and Signal Processing. Wiley-Interscience, 2000. [8] E. C. Ifeachor, Digital Signal Processing: A Practical Approach, second edition ed. Prentice-Hall, 2002. [9] A. Bateman and I. Paterson-Stephens, The DSP Handbook: Algorithms, Applications and Design Techniques. Prentice-Hall, 2002. [10] J. H. McClelland, R. W. Schafer, and M. A. Yoder, DSP First: A Multimedia Approach. Prentice-Hall, 1998. [11] P. S. Diniz, E. A. da Silva, and S. L. Netto, Digital Signal Processing: System Analysis and Design. Cambridge University Press, 2002.
208
BIBLIOGRAPHY
[12] S. Salivahanan, A. Vallavaraj, and C. Gnanapriya, Digital Signal Processing. Tata McGraw-Hill, 2000. [13] M. H. Hayes, Schaums Outline of Theory and Problems of Digital Signal Processing. McGraw-Hill, 1999. [14] T. K. Moon and W. C. Sterling, Mathematical Methods and Algorithms for Signal Processing. Prentice-Hall, 2000.
Index
ADPCM, 7 aliasing, 100, 102 arithmetic overow, 155, 167 autocorrelation, 56 beamforming, 15 bilinear transformation, 196 block processing, 4 bounded, 27 causal, 23 circular buffer, 153 circular convolution, 67 companding, 6, 122 -law, 123 A-law, 124 convolution sum, 37 correlation, 29 critical sampling, 101 DFT, see discrete Fourier Transform digital signals, 19 discrete Fourier transform, 17, 59 discrete-time Fourier series, 43, 44 discrete-time Fourier Transform, 17 discrete-time Fourier transform, 43, 49 convergence, 50 mean square, 51 uniform, 50 inverse, 56 properties, 53 frequency shift, 55 linearity, 54 periodicity, 53 symmetry, 53 time reversal, 55 time shift, 54 discrete-time signals, 19 discrete-time system frequency response, 134 minimum phase, 136 poles, 134 zeros, 134 Dithering, 124 down-sampling, 28 DTFT, see discrete-time Fourier transform echo cancellation, 12 energy density spectrum, 52 energy signal, 25 equalizer, 12 lter adaptive, 11 analog, 188 Butterworth, 188 Chebyshev, 188, 192 anti-aliasing, 107 anti-image postlter, 112 bandpass, 3 design, 182 design process, 171
210 equalizing, 12 high-pass, 3 design, 181 low-pass, 3 order, 41 realization, 137 canonical form, 142 cascade, 148 direct form, 139 FIR, 144 parallel, 149 transposed form, 143 recursive, 187 nite impulse reponse lter window method, 176 nite impulse response, 18, 40, 171 nite impulse response lter frequency sampling, 183 FIR, see nite impulse response xed-point, 155 oating point, 157 Fourier coefcients, 44 discrete-time, 45 frequency, 3 normalized, 47 frequency resolution computational, 92 frequency transformation, 199 frequency warping, 197 Gibbs phenomenon, 52 group delay, 135 ideal reconstruction, 109 IEEE 754, 157 IIR, see innite impulse response image compression, 10 image enhancement, 9
INDEX image processing, 9 image restoration, 10 impulse invariant method, 202 innite impulse response, 18, 40, 187 innite impulse response lter design, 188 Kotelnikov, 102 limit cycle, 165 linear convolution, 37, 55 properties, 38 linear phase, 136 linear predictive coding, 7 modulation, 29 multiply-and-accumulate, 139 noise cancellation, 11 noise-shaping, 121 Nyquist interval, 105 Nyquist rate, 101 over-sampling, 101, 119 ratio, 119 Parseval, 47 PCM, 7 phase delay, 135 power density spectrum, 48 power signal, 26 power spectral density, 118 pre-warping, 197 quantization, 97, 114 non-uniform, 122 quantization error, see quantization noise quantization noise, 116 quantization width, 115 quantizer
INDEX full-scale range, 115 quantizer resolution, 115 round-off, 155 round-off noise, 162 sample, 19 sample-by-sample processing, 4 sampling, 19, 97 frequency, 99 period, 98 rate, 99 uniform, 98 sampling frequency, 19 scaling, 27 sequence average power, 26, 46 bounded, 33 causal, 23 complex exponential, 21 energy, 25 nite length, 23 periodic, 24 random, 22 sinusoidal, 21 symmetry, 23 unit impulse, 20 unit step, 20 Shannon, 101 signal, 2 chirp, 93 discrete-time, 19 non-stationary, 93 reconstruction, 108 signal-to-quantization-noise ratio, 118 smart antennas, 15 software radio, 16 spectral analysis, 84
211 spectrogram, 93 spectrum, 3 leakage, 86 time-varying, 93 speech coding, 5 speech processing, 5 speech recognition, 9 speech synthesis, 9 staircase reconstructor, 111 superposition, 36 superposition principle, 30 system causal, 32 discrete-time, 30 linear, 30 linear time invariant, 31 linear time-invariant, 34 lossless, 34 shift invariant, 31 stability, 33 bounded-input bounded-output, 33 exponential, 33 truncation, 155 under-sampling, 101 uniform sampling theorem, 101, 108 unit delay, 28 up-sampling, 28 Whittaker, 102 window Blackman, 180 Hamming, 91, 176 main lobe, 91 Hanning, 180 Kaiser, 180 rectangular, 88, 176, 177
212 main lobe, 88 sidelobe, 88, 90 wordlength, 154 Z transform, 131 inverse, 133 properties, 132 region of convergence, 131 zero-padding, 64
INDEX

An Introduction To Digital Signal Processing: Edmund M-K. Lai

Uploaded by

Copyright:

Available Formats

An Introduction To Digital Signal Processing: Edmund M-K. Lai

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

An Introduction To Digital Signal Processing: Edmund M-K. Lai

Uploaded by

Copyright:

Available Formats

An Introduction to Digital Signal Processing

Edmund M-K. Lai 17th November 2003

Image and Video Processing . . . . . . . . . . . . . . . . . . . .

Image Restoration . . . . . . . . . . . . . . . . . . . . 10 Image Compression and Coding . . . . . . . . . . . . 10

iv 1.4.3.1 1.4.3.2 1.4.3.3 1.4.4 1.4.5 1.4.6 1.5 2

CONTENTS Noise Cancellation . . . . . . . . . . . . . . . . . . . . 11 Echo Cancellation . . . . . . . . . . . . . . . . . . . . 11 Channel Equalization . . . . . . . . . . . . . . . . . . 12

Control Applications . . . . . . . . . . . . . . . . . . . . . . . . 13 Sensor or Antenna Array Processing . . . . . . . . . . . . . . . 15

Digital Communication Transceivers . . . . . . . . . . . . . . . 16

Brief Overview of Chapters . . . . . . . . . . . . . . . . . . . . . . . . . 16 19

Discrete-time Signals and Systems 2.1

Types of Sequences . . . . . . . . . . . . . . . . . . . . . . . . . 22 Some Basic Operations on Sequences . . . . . . . . . . . . . . . 27

FIR and IIR Systems . . . . . . . . . . . . . . . . . . . . . . . . 40

CONTENTS 3 Discrete-time Fourier Analysis 3.1 3.2

3.3.4 3.3.5 3.4

Inverse DTFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 Autocorrelation Sequence . . . . . . . . . . . . . . . . . . . . . 56

4 Discrete Fourier Transform and Its Computation 4.1

vi 4.1.3 4.1.4 4.1.5 4.1.6 4.1.7 4.1.8 4.1.9 4.2

Sampling, Reconstruction and Quantization 5.1

Reconstruction of Analog Signals . . . . . . . . . . . . . . . . . . . . . 108 5.2.1 Ideal Reconstructor . . . . . . . . . . . . . . . . . . . . . . . . . 109

CONTENTS 5.2.2 5.2.3 5.3

vii Staircase Reconstructor . . . . . . . . . . . . . . . . . . . . . . . 111 Anti-image Postlters . . . . . . . . . . . . . . . . . . . . . . . 111

Dithering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 131

Digital Filters and Their Realizations 6.1

Finite Wordlength Effects . . . . . . . . . . . . . . . . . . . . . . . . . . 154

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169 171

FIR Filter Design 7.1 7.2

7.3 7.4 7.5 7.6 8

IIR Filter Design 8.1 8.2

CONTENTS 8.2.2 8.3

ix Chebyshev Filters . . . . . . . . . . . . . . . . . . . . . . . . . . 192

Impulse Invariant Method . . . . . . . . . . . . . . . . . . . . . . . . . . 202 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 206 207 209

Advantages of Processing Signals Digitally

Denition of Some Terms

Denition of Some Terms

1.3.1 Sample-by-sample Processing

1.4.1 Speech and Audio Processing

Step Size Adjustment

Figure 1.2: Block diagram of an adaptive differential pulse code modulator

Impulse Train For voiced speech

Vocal Tract Model Timevarying Filter Gain Synthesized Speech

Random Noise For unvoiced speech

Figure 1.3: A model of speech production

Application Areas 1.4.1.3 Speech Synthesis

1.4.2 Image and Video Processing

Image Compression and Coding

1.4.3 Adaptive Filtering

Primary signal from abdominal sensor

Adjustable Filter Filter coefficient adjustment Adaptive Algorithm

4 wire to 2 wire junction hybrid

Nearend signal and echo of farend signal

Adjustable Filter Filter coefficient adjustment Adaptive Algorithm

Cleaned nearend signal

Adjustable Filter (Equalizer)

Figure 1.6: An adaptive equaliser in training mode

Figure 1.7: A digital closed loop control system

1.4.5 Sensor or Antenna Array Processing

1.4.6 Digital Communication Transceivers

Brief Overview of Chapters

An overview of the remaining chapters in this book are as follows: