Digital Signal Processing/Printable version
This is the print version of Digital Signal Processing You won't see this message or any elements not part of the book's content when you print or preview this page. 
The current, editable version of this book is available in Wikibooks, the opencontent textbooks collection, at
https://en.wikibooks.org/wiki/Digital_Signal_Processing
Introduction
Digital Signal Processing (DSP) is the new technological revolution. edit
Who is This Book For? edit
What Will This Book Cover? edit
How is This Book Arranged edit
Software tool
Digital Signal Processing is a field of study that combines both mathematical theory and physical implementation. It makes no sense to consider a digital system without asking oneself how it will be implemented. In the design and analysis phase, some generalpurpose signal processing tools are available.
Matlab edit
MATLAB is an excellent (although expensive) tool for simulating systems, and for creating the evervaluable "proof of concept". This book will make several references to MATLAB, but don't get confused: This book will not teach how to program in MATLAB. If you would like to learn MATLAB, check out the book MATLAB Programming.
There are other free alternatives to MATLAB, with varying degrees of code compatibility.
Octave edit
GNU Octave is a free and Open Source alternative to MATLAB. Octave can be obtained from http://www.octave.org. It endeavors to be MATLABcompatible and largely it is (a lot of MATLAB code can be run using Octave, and vice versa), though some functionality is missing. For more information on Octave, see the Wikibook Octave Programming.
SciPy edit
SciPy is a Pythonbased set of libraries which allow to perform numeric calculations. As the preceding tools, it features a signal processing toolbox. Also, the python scripts can make use of matplotlib, a plotting library whose basic commands are very similar to MATLAB's.
Discrete Data
Continuous data is something that most people are familiar with. A quick analogy: when an analog data sensor (e.g., your eye) becomes active (you blink it open), it starts receiving input immediately (you see sunshine), it converts the input (optical rays) to a desired output (optic nerve signals), and sends the data off to its destination (your brain). It does this without hesitation, and continues doing so until the sensor turns off (you blink your eyes closed). The output is often called a data "stream"; once started, it might run forever, unless something tells it to stop. Now, instead of a physical sensor, if we're able to define our data mathematically in terms of a continuous function, we can calculate our data value at any point along the data stream. It's important to realize that this provides the possibility of an infinite (∞) number of data points, no matter how small the interval might be between the start and stop limits of the data stream.
This brings us to the related concept of Discrete data. Discrete data is noncontinuous, only existing at certain points along an input interval, and thereby giving us a finite number of data points to deal with. The analogy with the vision would be the illumination with a stroboscope. The scene viewed consists of a series of images. As all information between two images is lost, the frequency of the stroboscopic illumination should be high enough not to miss the movements of a moving object. This data can also be defined by a mathematical function, but one that is limited and can only be evaluated at the discrete points of input. These are called "discrete functions" to distinguish them from the continuous variety.
Discrete functions and data give us the advantage of being able to deal with a finite number of data points.
Sets and Series edit
Discrete data is displayed in sets such as:
X[n] = [1 2 3& 4 5 6]
We will be using the "&" symbol to denote the data item that occurs at point zero. Now, by filling in values for n, we can select different values from our series:
X[0] = 3 X[2] = 1 X[3] = 6
We can move the zero point anywhere in the set that we want. It is also important to note that we can pad a series with zeros on either side, so long as we keep track of the zeropoint:
X[n] = [0 0 0 0 1 2 3& 4 5 6] = [1 2 3& 4 5 6]
In fact, we assume that any point in our series without an explicit value is equal to zero. So if we have the same set:
X[n] = [1 2 3& 4 5 6]
We know that every value outside of our given range is zero:
X[100] = 0 X[100] = 0
Two Types of Discrete Data Sets edit
Data can be discrete in magnitude, discrete in time, or both. Here are some examples:
 Discrete in Time
 Discreteintime values only exist at specific points in time. If we try to take the value of a discreteintime data set at a time point where there is no data, our result will be zero. The image below shows an example waveform that is discrete in time, having gain values defined only at certain points along the time line. Notice that while it is discrete in time values, the gain value of each sample is not limited or quantized, and it may take any magnitude value.
 Discrete in Value
 Discreteinvalue series can only exist at certain values of magnitude, no matter when we look at them. For instance, we might say that a certain computer device can only handle integers, and no decimals. The image below shows a signal that is discrete in magnitude, but is not discrete in time. While the waveform indeed has a gain value at every point in time, the values are quantized to specific magnitudes, such as only integers. The steps between the allowed values produce the "staircase" effect in the image. If we try to take the value at a time point exactly when the magnitude transition occurs, a quantization rule must be in place to deal with it, setting the value to the appropriate previous or subsequent magnitude.
Stem Plots edit
Discrete data is frequently represented with a stem plot. Stem plots mark data points with dots, and draw a vertical line between the taxis (the horizontal time axis) and the dot:
F[n] = [5& 4 3 2 1]
About the Notation edit
The notation we use to denote the zero point of a discrete set was chosen arbitrarily. Textbooks on the subject will frequently use arrows, periods, or underscores to denote the zero position of a set. Here are some examples:
 v Y[n] = [1 2 3 4 5]
. Y[n] = [1 2 3 4 5]
_ Y[n] = [1 2 3 4 5]
All of these things are too tricky to write in wikibooks, so we have decided to use an ampersand (&) to denote the zeropoint. The ampersand is not used for any other purpose in this book, so hopefully we can avoid some confusion.
Sampling edit
Sampling is the process of converting continuous data into discrete data. The sampling process takes a snapshot of the value of a given input signal, rounds if necessary (for discreteinvalue systems), and outputs the discrete data. A common example of a sampler is an Analog to Digital Converter (ADC).
Let's say we have a function based on time (t). We will call this continuoustime function f(t):
Where u(t) is the unit step function. Now, if we want to sample this function, mathematically, we can plug in discrete values for t, and read the output. We will denote our sampled output as F[n]:
F[n] = 0 : n ≤ 0 F[1] = f(1) = 2 F[2] = f(2) = 4 F[100] = f(100) = 200
This means that our output series for F[n] is the following:
F[n] = [0& 2 4 6 8 10 12 ...]
Reconstruction edit
We digitize (sample) our signal, we do some magical digital signal processing on that signal, and then what do we do with the results? Frequently, we want to convert that digital result back into an analog signal. The problem is that the sampling process loses a lot of data. Specifically, all the data between the discrete data points is lost, and if the sampler rounds the value, then a small amount of error is built into the system that can never be recovered. A common example of a reconstructor is a Digital to Analog Converter (DAC).
Interpolation edit
When converting a digital signal into an analog signal, frequently a process called Interpolation is used to make the analog version a more likely representation of the signal. Interpolation essentially "connects the dots" between discrete data points, and the reconstructor then outputs an analog waveform with the dots connected. We will show the process below:
We start out with an uninterpolated stem plot:
We connect the points in our stem plot with straight lines (dotted lines):
We draw new points on that line, midway between our existing lines (dashed lines):
Discrete Operations
There are a number of different operations that can be performed on discrete data sets.
Arithmetic edit
Let's say that we have 2 data sets, A[n] and B[n]. We will define these two sets to have the following values:
A[n] = [W X& Y Z] B[n] = [I J& K L]
Arithmetic operations on discrete data sets are performed on an itembyitem basis. Here are some examples:
Addition edit
A[n] + B[n] = [(W+I) (X+J)& (Y+K) (Z+L)]
If the zero positions don't line up (like they conveniently do in our example), we need to manually line up the zero positions before we add.
Subtraction edit
A[n]  B[n] = [(WI) (XJ)& (YK) (ZL)]
Same as addition, only we subtract.
Multiplication edit
In this situation, we are using the "asterisk" (*) to denote multiplication. Realistically, we should use the " " symbol for multiplication. In later examples, we will use the asterisk for other operations.
A[n] * B[n] = [(W*I) (X*J)& (Y*K) (Z*L)]
Division edit
A[n] / B[n] = [(W/I) (X/J)& (Y/K) (Z/L)]
Even though we are using the square brackets for discrete data sets, they should NOT be confused with matrices. These are not matrices, and we do not perform matrix operations on these sets.
TimeShifting edit
If we timeshift a discrete data set, we essentially are moving the zerotime reference point of the set. The zero point represents "now", creates the starting point of our view of the data, and the point's location is typically established by a need of the processing we're involved in.
Let's say we have the data set F[n] with the values:
F[n] = [1 2 3& 4 5 ]
Then we can shift this set backwards by 1 data point as such:
F[n1] = [1 2& 3 4 5]
We can shift the data forward by 1 data point in a similar manner:
F[n+1] = [1 2 3 4& 5]
Discrete data values are time oriented. Values to the right of the zero point are called "future values", and values to the left are "past values". When the data set is populated on both sides of the zero reference, "future" and "past" are synthetic terms relative to the zero point "now", and don't refer to current physical time. In this context, the data values have already been collected and, as such, can only come from our past.
It is important to understand that a physicallyrealizable digital system that is still receiving data cannot perform operations on future values. It makes no sense to require a future value in order to make a calculation, because such systems are often getting input from a sensor in realtime, and future values simply don't exist.
Time Inversion edit
Let's say that we have the same data set, F[n]:
F[n] = [1 2 3& 4 5]
We can invert the data set as such:
F[n] = [5 4 3& 2 1]
We keep the zero point in the same place, and we flip all the data items around. in essence, we are taking a mirror image of the data set, and the zero point is the mirror.
Convolution edit
Convolution is a much easier operation in the discrete time domain than in the continuous time domain. Let's say we have two data sets, A[n] and B[n]:
A[n] = [1& 0 1 2] B[n] = [2& 2 2 1]
We will denote convolution with the asterisk (*) in this section; it isn't "multiplication" here. Our convolution function is shown like this:
Y[n] = A[n] * B[n]
and it specifies that we will store our convolution values in the set named Y[n].
Convolution is performed by following a series of steps involving both sets of data points. First, we time invert one of the data sets. It doesn't matter which one, so we can pick the easiest choice:
A[n] = [1& 0 1 2] B[n] = [1 2 2 2&]
Next, we line up the data vertically so that only one data item overlaps:
A[n] > [1& 0 1 2] B[n] > [1 2 2 2&]
Now, we are going to zeropad both sets, making them equal length, putting zeros in open positions as needed:
A[n] > [0 0 0 1& 0 1 2] B[n] > [1 2 2 2& 0 0 0]
Now, we will multiply the contents of each column, and add them together:
A[n] > [0 0 0 1& 0 1 2]
B[n] > [1 2 2 2& 0 0 0]
Y[m] =
This gives us the first data point: Y[m] = [2].
Next, we need to time shift B[n] forward in time by one point, and do the same process: multiply the columns, and add:
A[n] > [0 0 0 1& 0 1 2]
B[n1] > [1 2 2 2& 0 0 0]
Y[m+1] =
Repeat the "time shift, multiply, and add" steps until no more data points overlap.
A[n] > [0 0 0 1& 0 1 2]
B[n2] > [1 2 2 2& 0 0 0]
Y[m+2] =
A[n] > [0 0 0 1& 0 1 2]
B[n3] > [1 2 2 2& 0 0 0]
Y[m+3] =
A[n] > [0 0 0 1& 0 1 2]
B[n4] > [1 2 2 2& 0 0 0]
Y[m+4] =
A[n] > [0 0 0 1& 0 1 2]
B[n] > [1 2 2 2& 0 0 0]
Y[m+5] =
A[n] > [0 0 0 1& 0 1 2]
B[n+6] > [1 2 2 2& 0 0 0]
Y[m+6] =
Now we have our full set of data points, Y[m]:
Y[m] = [2 2 4 7 6 5 2]
We have our values, but where do we put the zero point? It turns out that the zero point of the result occurs where the zero points of the two operands overlapped in our shift calculations. The result set becomes:
Y[n] = [2& 2 4 7 6 5 2]
It is important to note that the length of the result set is the sum of the lengths of the unpadded operands, minus 1. Or, if you prefer, the length of the result set is the same as either of the zeropadded sets (since they are equal length).
Discrete Operations in Matlab/Octave edit
These discrete operations can all be performed in Matlab, using special operators.
 Division and Multiplication
 Matlab is matrixbased, and therefore the normal multiplication and division operations will be the matrix operations, which are not what we want to do in DSP. To multiply items on a peritem basis, we need to prepend the operators with a period. For instance:
Y = X .* H %X times H Y = X ./ H %X divided by H
If we forget the period, Matlab will attempt a matrix operation, and will alert you that the dimensions of the matrixes are incompatible with the operator.
 Convolution
 The convolution operation can be performed with the conv command in Matlab. For instance, if we wanted to compute the convolution of X[n] and H[n], we would use the following Matlab code:
Y = conv(X, H);
or
Y = conv(H, X);
Difference Calculus edit
When working with discrete sets, you might wonder exactly how we would perform calculus in the discrete time domain. In fact, we should wonder if calculus as we know it is possible at all! It turns out that in the discrete time domain, we can use several techniques to approximate the calculus operations of differentiation and integration. We call these techniques Difference Calculus.
Derivatives edit
What is differentiation exactly? In the continuous time domain, the derivative of a function is the slope of the function at any given point in time. To find the derivative, then, of a discrete time signal, we need to find the slope of the discrete data set.
The data points in the discrete time domain can be treated as geometrical points, and we know that any two points define a unique line. We can then use the algebraic equation to solve for the slope, m, of the line:
Now, in the discrete time domain, f(t) is sampled and replaced by F[n]. We also know that in discrete time, the difference in time between any two points is exactly 1 unit of time! With this mind, we can plug these values into our above equation:
or simply:
We can then timeshift the entire equation one point to the left, so that our equation doesn't require any future values:
The derivative can be found by subtracting timeshifted (delayed) versions of the function from itself.
Difference Calculus Equations edit
Difference calculus equations are arithmetic equations, with a few key components. Let's take an example:
We can see in this equation that X[n1] has a coefficient of 2, and X[n2] has a coefficient of 5. In difference Calculus, the coefficients are known as "tap weights". We will see why they are called tap weights in later chapters about digital filters.
Digital Systems
Systems Introduction edit
Digital systems can be conceptually very difficult to understand at first. Let's start out with a block diagram:
We have a digital system, h[n], that is going to alter the input (x[n]) in some way, and produce an output (y[n]). At each integer value for discrete time, n, we will feed 1 value of x[n] into the machine, turn the crank, and create 1 output value for y[n].
Basic Example edit
Let's say that h[n] is a multiplier circuit with a value of 5. Every input value is multiplied by 5, and that is the output. In other words, we can show our difference calculus equation as such:
now, for each successive value for n, we can calculate our output value, y[n]. As an example, using the above difference equation, we can feed in an experimental input:
x[n] = [1 0 1 2]
And by multiplying each data item by 5, we get the following result:
y[n] = [5 0 5 10]
Properties of Digital systems edit
A loose definition of a causal system is a system in which the output does not change before the input. All real systems are causal, and it is impossible to create a system that is not causal.
Circuit Symbols edit
 wires
 pickoff nodes
 adders
 multipliers
Physically Realizeable Systems edit
There is a distinction to be made between systems which are possible mathematically, and systems that can be physically implemented in digital hardware. Some systems might look good on paper, but they are difficult, if not impossible, to actually build. One criterion we have already seen is that a physically realizable digital system must not rely on future values, only on past values. We will discuss other criteria in future chapters, as the concepts involved become more familiar.
Impulse Response
Let's say that we have the following block diagram:
h[n] is known as the 'Impulse Response of the digital system. We will talk about impulse responses here in this chapter.
Impulse Function edit
In digital circuits, we use a variant of the continuoustime delta function. This delta function (δ[n]) has a value of 1 at n = 0, and a value of 0 everywhere else.
δ[n] = [... 0 0 0 0 1& 0 0 0 0 ...]
If we set x[n] = δ[n], then the output y[n] will be the response of the system to an impulse, or simply, the "Impulse Response".
We can timeshift the impulse function as such:
δ[n1] = [... 0 0 0 0& 1 0 0 0 0 ...] δ[n+1] = [... 0 0 0 0 1 0& 0 0 0 ...]
We can add timeshifted values of the impulse function to create an Impulse Train:
y[n] = δ[n] + δ[n1] + δ[n2] + δ[n4] y[n] = [1& 1 1 0 1]
Impulse Trains edit
Impulse Response edit
If we have a difference equation relating y[n] to x[n], we can find the impulse response difference equation by replacing every y with an h, and every x with a δ:
y[n] = 2x[n] + 3x[n1] + 4x[n2] h[n] = 2δ[n] + 3δ[n1] + 4δ[n2]
And by plugging in successive values for n, we can calculate the impulse response to be:
h[n] = [2& 3 4]
Output edit
Now, let's say that we have a given impulse response, h[n], and we have a given input x[n] as such:
x[n] = [1& 0 1 2] h[n] = [2& 2 2 1]
We can calculate the output, y[n], as the convolution of those 2 quantities:
x[n] > 1& 0 1 2 h[n] > 1 2 2 2& h[n3] > 1 2 2 2& > y[m3] = 2& h[n2] > 1 2 2 2& > y[m2] = 2 h[n1] > 1 2 2 2& > y[m1] = 4 h[n] > 1 2 2 2& > y[m] = 7 h[n+1] > 1 2 2 2& > y[m+1] = 6 h[n+2] > 1 2 2 2& > y[m+2] = 5 h[n+3] > 1 2 2 2& > y[m+3] = 2 y[n] = [2& 2 4 7 6 5 2]
Sampling and Reconstruction
digital pro{sampling} edit
Sampling is the process of recording the values of a signal at given points in time. For A/D converters, these points in time are equidistant. The number of samples taken during one second is called the sample rate. Keep in mind that these samples are still analogue values. The mathematic description of the ideal sampling is the multiplication of the signal with a sequence of direct pulses. In real A/D converters the sampling is carried out by a sampleandhold buffer. The sampleandhold buffer splits the sample period in a sample time and a hold time. In case of a voltage being sampled, a capacitor is switched to the input line during the sample time. During the hold time it is detached from the line and keeps its voltage.
Quantization edit
Quantization is the process of representing the analog voltage from the sampleandhold circuit by a fixed number of bits. The input analog voltage (or current) is compared to a set of predefined voltage (or current) levels. Each of the levels is represented by a unique binary number, and the binary number that corresponds to the level that is closest to the analog voltage is chosen to represent that sample. This process rounds the analog voltage to the nearest level, which means that the digital representation is an approximation to the analog voltage. There are a few methods for quantizing samples. The most commonly used ones are the dual slope method and the successive approximation.
Dual Slope edit
Successive Approximation edit
Reconstruction edit
Reconstruction is the process of creating an analog voltage (or current) from samples. A digitaltoanalog converter takes a series of binary numbers and recreates the voltage (or current) levels that corresponds to that binary number. Then this signal is filtered by a lowpass filter. This process is analogous to interpolating between points on a graph, but it can be shown that under certain conditions the original analog signal can be reconstructed exactly from its samples. Unfortunately, the conditions for exact reconstruction cannot be achieved in practice, and so in practice the reconstruction is an approximation to the original analog signal.
Aliasing edit
Aliasing is a common problem in digital media processing applications. Many readers have heard of "antialiasing" features in highquality video cards. This page will explain what Aliasing is, and how it can be avoided.
Aliasing is an effect of violating the NyquistShannon sampling theory. During sampling the base band spectrum of the sampled signal is mirrored to every multifold of the sampling frequency. These mirrored spectra are called alias. If the signal spectrum reaches farther than half the sampling frequency base band spectrum and aliases touch each other and the base band spectrum gets superimposed by the first alias spectrum. The easiest way to prevent aliasing is the application of a steep sloped lowpass filter with half the sampling frequency before the conversion. Aliasing can be avoided by keeping Fs>2Fmax.
Nyquist Sampling Rate edit
The Nyquist Sampling Rate is the lowest sampling rate that can be used without having aliasing. The sampling rate for an analog signal must be at least two times the bandwidth of the signal. So, for example, an audio signal with a bandwidth of 20 kHz must be sampled at least at 40 kHz to avoid aliasing. In audio CD's, the sampling rate is 44.1 kHz, which is about 10% higher than the Nyquist Sampling Rate to allow cheaper reconstruction filters to be used. The Nyquist Sampling Rate is the lowest sampling rate that can be used without having aliasing.
AntiAliasing edit
The sampling rate for an analog signal must be at least two times as high as the highest frequency in the analog signal in order to avoid aliasing. Conversely, for a fixed sampling rate, the highest frequency in the analog signal can be no higher than a half of the sampling rate. Any part of the signal or noise that is higher than a half of the sampling rate will cause aliasing. In order to avoid this problem, the analog signal is usually filtered by a low pass filter prior to being sampled, and this filter is called an antialiasing filter. Sometimes the reconstruction filter after a digitaltoanalog converter is also called an antialiasing filter.
Converters edit
As a matter of professional interest, we will use this page to briefly discuss AnalogtoDigital (A2D) and DigitaltoAnalog (D2A) converters, and how they are related to the field of digital signal processing. Strictly speaking, this page is not important to the core field of DSP, so it can be skipped at the reader's leisure.
A2D converters are the bread and butter of DSP systems. A2D converters change analog electrical data into digital signals through a process called sampling.
D2A converters attempt to create an analog output waveform from a digital signal. However, it is nearly impossible to create a smooth output curve in this manner, so the output waveform is going to have a wide bandwidth, and will have jagged edges. Some techniques, such as interpolation can be used to improve these output waveforms.
ContinuousTime Fourier Transform
The ContinuousTime Fourier Transform (CTFT) is the version of the fourier transform that is most common, and is the only fourier transform so far discussed in EE wikibooks such as Signals and Systems, or Communication Systems. This transform is mentioned here as a stepping stone for further discussions of the DiscreteTime Fourier Transform (DTFT), and the Discrete Fourier Transform (DFT). The CTFT itself is not useful in digital signal processing applications, so it will not appear much in the rest of this book
CTFT Definition edit
The CTFT is defined as such:
[CTFT]
CTFT Use edit
Frequency Domain edit
Convolution Theorem edit
Multiplication in the time domain becomes convolution in the frequency domain. Convolution in the time domain becomes multiplication in the frequency domain. This is an example of the property of duality. The convolution theorem is a very important theorem, and every transform has a version of it, so it is a good idea to become familiar with it now (if you aren't previously familiar with it).
Further reading edit
For more information about the Fourier Transform, see Signals and Systems.
DiscreteTime Fourier Transform
The DiscreteTime Fourier Transform (DTFT) is the cornerstone of all DSP, because it tells us that from a discrete set of samples of a continuous function, we can create a periodic summation of that function's Fourier transform. At the very least, we can recreate an approximation of the actual transform and its inverse, the original continuous function. And under certain idealized conditions, we can recreate the original function with no distortion at all. That famous theorem is called the NyquistShannon sampling theorem.
DTFT edit
[DTFT]
The resulting function, is a continuous function that is interesting for analysis. It can be used in programs, such as Matlab, to design filters and obtain the corresponding timedomain filter values.
DTFT Convolution Theorem edit
Like the CTFT, the DTFT has a convolution theorem associated with it. However, since the DTFT results in discretefrequency values, the convolution theorem needs to be modified as such:
 DTFT Convolution Theorem
 Multiplication in the continuous time domain becomes discrete convolution in the discrete frequency domain. Convolution in the continuous time domain becomes multiplication in the discrete frequency domain.
Energy edit
It is sometimes helpful to calculate the amount of energy that exists in a certain set. These calculations are based off the assumption that the different values in a set are voltage values, however this doesn't necessarily need to be the case to employ these operations.
We can show that the energy of a given set can be given by the following equation:
Energy in Frequency edit
Likewise, we can make a formula that represents the power in the continuousfrequency output of the DTFT:
Parseval's Theorem edit
Parseval's theorem states that the energy amounts found in the time domain must be equal to the energy amounts found in the frequency domain:
Power Density Spectrum edit
We can define the power density spectrum of the continuoustime frequency output of the DTFT as follows:
The area under the power density spectrum curve is the total energy of the signal.
Discrete Fourier Transform
As the name implies, the Discrete Fourier Transform (DFT) is purely discrete: discretetime data sets are converted into a discretefrequency representation. This is in contrast to the DTFT that uses discrete time, but converts to continuous frequency. Since the resulting frequency information is discrete in nature, it is very common for computers to use DFT(Discrete fourier Transform) calculations when frequency information is needed.
Using a series of mathematical tricks and generalizations, there is an algorithm for computing the DFT that is very fast on modern computers. This algorithm is known as the Fast Fourier Transform (FFT), and produces the same results as the normal DFT, in a fraction of the computational time as ordinary DFT calculations. The FFT will be discussed in detail in later chapters.
DFT edit
The Discrete Fourier Transform is a numerical variant of the Fourier Transform. Specifically, given a vector of n input amplitudes such as {f_{0}, f_{1}, f_{2}, ... , f_{n2}, f_{n1}}, the Discrete Fourier Transform yields a set of n frequency magnitudes.
The DFT is defined as such:
here, k is used to denote the frequency domain ordinal, and n is used to represent the timedomain ordinal. The big "N" is the length of the sequence to be transformed.
The Inverse DFT (IDFT) is given by the following equation:
Where is defined as:
Again, "N" is the length of the transformed sequence.
Matrix Calculations edit
The DFT can be calculated much more easily using a matrix equation:
The little "x" is the time domain sequence arranged as a Nx1 vertical vector. The big "X" is the resulting frequency information, that will be arranged as a vertical vector (as per the rules of matrix multiplication). The "D_{N}" term is a matrix defined as such:
In a similar way, the IDFT can be calculated using matrices as follows:
And we will define as the following:
Visualization of the Discrete Fourier Transform edit
It may be easily seen that the term
represents a unit vector in the complex plane, for any value of j and k. The angle of the vector is initially 0 radians (along the real axis) for or . As j and k increase, the angle is increased in units of 1/n^{th} of a circle. The total angle therefore becomes radians.
To understand the meaning of the Discrete Fourier Transform, it becomes effective to write the transform in matrix form, depicting the complex terms pictorally.
For example, with , the fourier transform can be written:
The top row (or the lefthand column) changes not at all; therefore, F0 is the resultant when the amplitudes augment each other. The second row changes by 1/n^{th} of a circle, going from left to right. The third row changes by 2/n^{th}s of a circle; the fourth row changes by 3/n^{th}s of a circle; the fifth row, here, changes by 4/n^{th}s = 4/8 = 1/2 a circle. Therefore, F4 is the resultant when the even terms are set at odds with the odd terms.
We therefore see from the square matrix that the DFT is the resultant of all 2D vectoral combinations of the input amplitudes.
Relation to DTFT edit
The DFT and the DTFT are related to each other in a very simple manner. If we take the DTFT of a given time sequence, x[n], we will get a continuousfrequency function called . If we sample at N equallyspaced locations, the result will be the DFT, X[k]
Circular Time Shifting edit
Circular Time Shifting is very similar to regular, linear time shifting, except that as the items are shifted past a certain point, they are looped around to the other end of the sequence. This subject may seem like a bit of a tangent, but the importance of this topic will become apparent when we discuss the Circular Convolution operation in the next chapter.
Let's say we have a given sequence, x[n], as such:
x[n] = [1& 2 3 4]
As we know, we can timeshift this sequence to the left or the right by subtracting or adding numbers to the argument:
x[n+1] = [1 2& 3 4 0]
x[n1] = [0& 1 2 3 4]
Here, we've padded both sequences with zeros to make it obvious where the next item is coming in from. Now, let's say that when we shift a set, instead of padding with zeros, we loop the first number around to fill the hole. We will use the notation x[<n+A>] to denote an Apoint circular shift:
x[<n+1>] = [2& 3 4 1]
Notice how the length of the sequence started as 4, and ended as 4, and we didn't need to pad with zeros. Also notice that the 1 "fell off" the left side of the set, and "reappeared" on the right side. It's as if the set is a circle, and objects that rotate off one direction will come back from the other direction. Likewise, let us rotate in the other direction:
x[<n1>] = [4& 1 2 3]
If we rotate around by N points, where N is the length of the set, we get an important identity:
x[<n+N>] = x[<nN>] = x[n]
Time Inversion edit
Now, let's look at a more complicated case: mixing circular shifting with time inversion. Let's define our example set x[n]:
x[n] = [1& 2 3 4]
we can define our timeinverted set as such:
x[n] = [4 3 2 1&]
notice that the "4" went from n=3 to n=3, but the "1" stayed at the zero point. Now, let's try to combine this with our circular shift operation, to produce a result x[<n>]. We know that the "1" (in the zero position) has not moved in either direction, nor have we shifted to the left or right, because we aren't adding or subtracting anything to our argument (n). So we can put 1 back in the same place as it was. All the rest of the numbers (2 3 and 4) all shifted to the right of the sequence, and therefore we will loop them around to the other side:
x[<n>] = [1& 4 3 2]
Here is a general rule for circular timeinversion:
Circular Convolution edit
Circular Convolution is an important operation to learn, because it plays an important role in using the DFT.
Let's say we have 2 discrete sequences both of length N, x[n] and h[n]. We can define a circular convolution operation as such:
notice how we are using a circular timeshifting operation, instead of the linear timeshift used in regular, linear convolution.
Circular Convolution Theorem edit
The DFT has certain properties that make it incompatible with the regular convolution theorem. However, the DFT does have a similar relation with the Circular Convolution operation that is very helpful. We call this relation the Circular Convolution Theorem, and we state it as such:
 Circular Convolution Theorem
 Multiplication in the discretetime domain becomes circular convolution in the discretefrequency domain. Circular convolution in the discretetime domain becomes multiplication in the discretefrequency domain.
Relation to Linear Convolution edit
Circular Convolution is related to linear convolution, and we can use the circular convolution operation to compute the linear convolution result. Let's say we have 2 sets, x[n] and h[n], both of length N:
 pad both sets with N1 zeros on the right side.
 perform the circular convolution on the new sets of length N+(N1)
 The result is the linear convolution result of the original two sets.
Fast Fourier Transform (FFT) edit
The Fast Fourier Transform refers to algorithms that compute the DFT in a numerically efficient manner.Algorithms like:decimation in time and decimation in frequency. This is a algorithm for computing the DFT that is very fast on modern computers.
Z Transform
The Z Transform has a strong relationship to the DTFT, and is incredibly useful in transforming, analyzing, and manipulating discrete calculus equations. The Z transform is named such because the letter 'z' (a lowercase Z) is used as the transformation variable.
z Transform Definition edit
For a given sequence x[n], we can define the ztransform X(z) as such:
[Z Transform]
 it is important to note that z is a continuous complex variable defined as such:
 where is the imaginary unit.
There can be several sequences which will generate the same ztransform with the different functions being differentiated by the convergence region of for which the summation in the ztransform will converge. These convergence regions are annuli centered at the orgin. In a given convergence region, only one will converge to a given .
Example 1 edit
 for
Example 2 edit
 for
Note that both examples have the same function as their ztransforms but with different convergence regions needed for the respective infinite summations in their ztransforms to converge. Many textbooks on ztransforms are only concerned with socalled righthanded functions which is to say functions for which for all less than some initial start point ; that is, for all . So long as the function grows at most exponentially after the start point, the ztransform of these socalled righthanded functions will converge in an open annulus going to infinity, for some positive real .
It is important to note that the ztransform rarely needs to be computed manually, because many common results have already been tabulated extensively in tables, and control system software includes it (MatLab,Octave,SciLab).
The ztransform is actually a special case of the socalled Laurent series, which is a special case of the commonly used Taylor series.
The Inverse ZTransform edit
The inverse ztransform can be defined as such:
[Inverse Z Transform]
Where C is a closedcontour that lies inside the unit circle on the zplane, and encircles the point z = {0, 0}.
The inverse ztransform is mathematically very complicated, but luckily—like the ztransform itself—the results are extensively tabulated in tables.
Equivalence to DTFT edit
If we substitute , where is the frequency in radians per second, into the Ztransform, we get
which is equivalent to the definition of the DiscreteTime Fourier Transform. In other words, to convert from the Ztransform to the DTFT, we need to evaluate the Ztransform around the unit circle.
Properties edit
Since the ztransform is equivalent to the DTFT, the ztransform has many of the same properties. Specifically, the ztransform has the property of duality, and it also has a version of the convolution theorem (discussed later).
The ztransform is a linear operator.
Convolution Theorem edit
Since the Ztransform is equivalent to the DTFT, it also has a convolution theorem that is worth stating explicitly:
 Convolution Theorem
 Multiplication in the discretetime domain becomes convolution in the zdomain. Multiplication in the zdomain becomes convolution in the discretetime domain.
Y(s)=X(s).H(s)
ZPlane edit
Since the variable z is a continuous, complex variable, we can map the z variable to a complex plane as such:
Transfer Function edit
Let's say we have a system with an input/output relationship defined as such:
Y(z) = H(z)X(z)
We can define the transfer function of the system as being the term H(z). If we have a basic transfer function, we can break it down into parts:
Where H(z) is the transfer function, N(z) is the numerator of H(z) and D(z) is the denominator of H(z). If we set N(z)=0, the solutions to that equation are called the zeros of the transfer function. If we set D(z)=0, the solutions to that equation are called the poles of the transfer function.
The poles of the transfer function amplify the frequency response while the zero's attenuate it. This is important because when you design a filter you can place poles and zero's on the unit circle and quickly evaluate your filters frequency response.
Example edit
Here is an example:
So by dividing through by X(z), we can show that the transfer function is defined as such:
We can also find the D(z) and N(z) equations as such:
And from those equations, we can find the poles and zeros:
 Zeros
 z → 0
 Poles
 z → 1/2
Stability edit
It can be shown that for any causal system with a transfer function H(z), all the poles of H(z) must lie within the unitcircle on the zplane for the system to be stable. Zeros of the transfer function may lie inside or outside the circle. See Control Systems/Jurys Test.
Gain edit
Gain is the factor by which the output magnitude is different from the input magnitude. If the input magnitude is the same as the output magnitude at a given frequency, the filter is said to have "unity gain".
Properties edit
Here is a listing of the most common properties of the Z transform.
Time domain  Zdomain  ROC  

Notation  ROC:  
Linearity  At least the intersection of ROC_{1} and ROC_{2}  
Time shifting  ROC, except if and if  
Scaling in the zdomain  
Time reversal  
Conjugation  ROC  
Real part  ROC  
Imaginary part  ROC  
Differentiation  ROC  
Convolution  At least the intersection of ROC_{1} and ROC_{2}  
Correlation  At least the intersection of ROC of X_{1}(z) and X_{2}( )  
Multiplication  At least  
Parseval's relation 
 Initial value theorem
 , If causal
 Final value theorem
 , Only if poles of are inside unit circle
Phase Shift edit
Further reading edit
Hilbert Transform
Definition edit
The Hilbert transform is used to generate a complex signal from a real signal. The Hilbert transform is characterized by the impulse response:
The Hilbert Transform of a function x(t) is the convolution of x(t) with the function h(t), above. The Hilbert Transform is defined as such:
[Hilbert Transform]
We use the notation to denote the Hilbert transformation of x(t).
Fast Fourier Transform (FFT) Algorithm
The Fourier Transform (FFT) is an algorithm used to calculate a Discrete Fourier Transform (DFT) using a smaller number of multiplications.
Cooley–Tukey algorithm edit
The Cooley–Tukey algorithm is the most common fast Fourier transform (FFT) algorithm. It breaks down a larger DFT into smaller DFTs. The best known use of the Cooley–Tukey algorithm is to divide a N point transform into two N/2 point transforms, and is therefore limited to poweroftwo sizes.
The algorithm can be implemented in two different ways, using the socalled:
 Decimation In Time (DIT)
 Decimation In Frequency (DIF)
Circuit diagram edit
The following picture represents a 16point DIF FFT circuit diagram:
In the picture, the coefficients are given by:
where N is the number of points of the FFT. They are evenly spaced on the lower half of the unit circle and .
Test script edit
The following python script allows to test the algorithm illustrated above:
#!/usr/bin/python3
import math
import cmath
# 
# Parameters
#
coefficientBitNb = 8 # 0 for working with reals
# 
# Functions
#
# ..............................................................................
def bit_reversed(n, bit_nb):
binary_number = bin(n)
reverse_binary_number = binary_number[1:1:1]
reverse_binary_number = reverse_binary_number + \
(bit_nb  len(reverse_binary_number))*'0'
reverse_number = int(reverse_binary_number, bit_nb1)
return(reverse_number)
# ..............................................................................
def fft(x):
j = complex(0, 1)
# sizes
point_nb = len(x)
stage_nb = int(math.log2(point_nb))
# vectors
stored = x.copy()
for index in range(point_nb):
stored[index] = complex(stored[index], 0)
calculated = [complex(0, 0)] * point_nb
# coefficients
coefficients = [complex(0, 0)] * (point_nb//2)
for index in range(len(coefficients)):
coefficients[index] = cmath.exp(j * index/point_nb * 2*cmath.pi)
coefficients[index] = coefficients[index] * 2**coefficientBitNb
# loop
for stage_index in range(stage_nb):
# print([stored[i].real for i in range(point_nb)])
index_offset = 2**(stage_nbstage_index1)
# butterfly additions
for vector_index in range(point_nb):
isEven = (vector_index // index_offset) % 2 == 0
if isEven:
operand_a = stored[vector_index]
operand_b = stored[vector_index + index_offset]
else:
operand_a =  stored[vector_index]
operand_b = stored[vector_index  index_offset]
calculated[vector_index] = operand_a + operand_b
# coefficient multiplications
for vector_index in range(point_nb):
isEven = (vector_index // index_offset) % 2 == 0
if not isEven:
coefficient_index = (vector_index % index_offset) \
* (stage_index+1)
# print( \
# "(%d, %d) > %d" \
# % (stage_index, vector_index, coefficient_index) \
# )
calculated[vector_index] = \
coefficients[coefficient_index] * calculated[vector_index] \
/ 2**coefficientBitNb
if coefficientBitNb > 0:
calculated[vector_index] = complex( \
math.floor(calculated[vector_index].real), \
math.floor(calculated[vector_index].imag) \
)
# storage
stored = calculated.copy()
# reorder results
for index in range(point_nb):
calculated[bit_reversed(index, stage_nb)] = stored[index]
return calculated
# 
# Main program
#
source = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]
transformed = fft(source)
amplitude = [0.0] * len(source)
amplitude_print = ''
for index in range(len(transformed)):
amplitude[index] = abs(transformed[index])
amplitude_print = amplitude_print + "%5.3f " % amplitude[index]
print()
print(amplitude_print)
It has a special parameter, coefficientBitNb
, which allows to determine the calculation results
for a circuit using only integer numbers.
Digital Filters
Digital Filters can be very complicated devices, but they must be able to map to the difference equations of the filter design. This means that since difference equations only have a limited number of operations available (addition and multiplication), digital filters only have limited operations that they need to handle as well. There are only a handful of basic components to a digital filter, although these few components can be arranged in complex ways to make complicated filters.
Digital Filter Types edit
There are two types of filters in the digital realm: Finite Impulse Response (FIR) filters and Infinite Impulse Response (IIR) filters. They are very different in essence.
FIR Filters edit
FIR filters are specific to sampled systems. There is no equivalent in continuoustime systems. Therefore, only very specific analog filters are capable of implementing an FIR filter.
If we define the discrete time impulse function as
the response of an FIR filter to δ[n], denoted as h[n], will decay to zero after a finite number of samples. The transfer function of an FIR filter contains only zeros and either no poles or poles only at the origin.
An FIR filter with symmetric coefficients is guaranteed to provide a linear phase response, which can be very important in some applications. However, FIR filters suffer from low efficiency, and creating an FIR to meet a given spec requires much more hardware than an equivalent IIR filter.
IIR Filters (infinite impulse response filter) edit
IIR filters are typically designed basing on continuoustime transfer functions. IIR filters differ from FIR filters because they always contain feedback elements in the circuit, which can make the transfer functions more complicated to work with.
The transfer function of an IIR filter contains both poles and zeros. Its impulse response never decays to zero (though it may get so close to zero that the response cannot be represented with the number of bits available in the system).
IIR filters provide extraordinary benefits in terms of computing: IIR filters are more than an order of magnitude more efficient than an equivalent FIR filter. Even though FIR is easier to design, IIR will do the same work with fewer components, and fewer components translate directly to less money.
Filtering a timeseries in Octave/Sciplot edit
First create some points for a time series. In this case we'll create one second of random data sampled at 44 kHz.
sampling_t = 1/44000; t = 0:sampling_t:1; x = rand(size(t));
Have a look at the time series
plot(t,x);
Have a look at its spectrum (it is mostly uniform, what we would expect from noise)
specgram(x)
Now we'll use a builtin function to create a third order Butterworth lowpass filter with cutoff frequency pi*0.1 radians
[b,a] = butter(3,0.1)
Now filter the time series.
y = filter(b,a,x);
Have a look at the first 100 points of the filtered data.
hold on plot(y(1:100)) plot(x(1:100)) hold off
Check its spectrogram
specgram(y)
Now look at the magnitude of the FFT
plot(log(abs(fft(y))))
Filter response Types edit
Highpass and Lowpass edit
HighPass and LowPass filters are the simplest forms of digital filters, and they are relatively easy to design to specifications. This page will discuss highpass and lowpass transfer functions, and the implementations of each as FIR and IIR designs.
Bandpass edit
Bandpass Filters are like a combination of a highpass and a lowpass filter. Only specific bands are allowed to pass through the filter. Frequencies that are too high or too low will be rejected by the filter.
Stopband edit
Notch filters are the complement of Bandpass filters in that they only stop a certain narrow band of frequency information, and allow all other data to pass without problem.
Notch edit
The direct complement of a bandpass filter is called a bandstop filter. A notch filter is essentially a very narrow bandstop filter.
Comb Filters edit
Comb Filters, as their name implies, look like a hair comb. They have many "teeth", which in essence are notches in the transfer function where information is removed. These notches are spaced evenly across the spectrum, so they are only useful for removing noise that appear at regular frequency intervals.
Allpass edit
The Allpass filter is a filter that has a unity magnitude response for all frequencies. It does not affect the gain response of a system. The Allpass filter does affect the phase response of the system. For example, if an IIR filter is designed to meet a prescribed magnitude response often the phase response of the system will become nonlinear. To correct for this nonlinearity an Allpass filter is cascaded with the IIR filter so that the overall response(IIR+Allpass) has a constant group delay.
Canonic and NonCanonic edit
Canonic filters are filters where the order of the transfer function matches the number of delay units in the filter. Conversely, Noncanonic filters are when the filter has more delay units than the order of the transfer function. In general, FIR filters are canonic. Some IIR filter implementations are noncanonic.
Here is an example of a canonical, 2nd order filter:
Here is an example of a noncanonical 2nd order filter:
Notice that in the canonical filter, the system order (here: 2) equals the number of delay units in the filter (2). In the noncanonical version, the system order is not equal to the number of delays (4 delay units). Both filters perform the same task. The canonical version results in smaller digital hardware, because it uses fewer delay units. However, this drawback nearly disappears if several second order IIRs are cascaded, as they can share delay elements. In this case, the only "extra" delay elements are those on the input side of the first section.
Notations edit
There are some basic notations that we will be using:
 ω_{p} Passband cut off frequency
 ω_{s} Stopband cut off frequency
 δ_{p} Passband ripple peak magnitude
 δ_{s} Stopband ripple peak magnitude
FIR Filter Design
Filter design edit
The design procedure most frequently starts from the desired transfer function amplitude. The inverse Laplace transform provides the impulse response in the form of sampled values.
The length of the impulse response is also called the filter order.
Ideal lowpass filter edit
The ideal (brickwall or sinc filter) lowpass filter has an impulse response in the form of a sinc function:
This function is of infinite length. As such it can not be implemented in practice. Yet it can be approximated via truncation. The corresponding transfer function ripples on the sides of the cutoff frequency transition step. This is known as the Gibbs phenomenon.
Window design method edit
In order to smoothout the transfer function ripples close to the cutoff frequency, one can replace the brickwall shape by a continuous function such as the Hamming, Hanning, Blackman and further shapes. These functions are also called window functions and are also used for smoothing a set of samples before processing it.
The following code illustrates the design of a Hamming window FIR:
#!/usr/bin/python3
import math
import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt
# 
# Constants
#
# filter
signal_bit_nb = 16
coefficient_bit_nb = 16
sampling_rate = 48E3
cutoff_frequency = 5E3
filter_order = 100
# time signal
input_signal_duration = 10E3
frequency_1 = 1E3
amplitude_1 = 0.5
frequency_2 = 8E3
amplitude_2 = 0.25
# display
plot_time_signals = False
plot_transfer_function = False
plot_zeros = True
input_input_signal_color = 'blue'
filtered_input_signal_color = 'red'
transfer_function_amplitude_color = 'blue'
transfer_function_phase_color = 'red'
locus_axes_color = 'deepskyblue'
locus_zeroes_color = 'blue'
locus_cutoff_frequency_color = 'deepskyblue'
#
# Filter design
#
Nyquist_rate = sampling_rate / 2
sampling_period = 1/sampling_rate
coefficient_amplitude = 2**(coefficient_bit_nb1)
FIR_coefficients = sig.firwin(filter_order, cutoff_frequency/Nyquist_rate)
FIR_coefficients = np.round(coefficient_amplitude * FIR_coefficients) \
/ coefficient_amplitude
transfer_function = sig.TransferFunction(
FIR_coefficients, [1],
dt=sampling_period
)
# 
# Time filtering
#
sample_nb = round(input_signal_duration * sampling_rate)
# time signal
time = np.arange(sample_nb) / sampling_rate
# input signal
input_signal = amplitude_1 * np.sin(2*np.pi*frequency_1*time) \
+ amplitude_2*np.sin(2*np.pi*frequency_2*time)
input_signal = np.round(2**(signal_bit_nb1) * input_signal)
# filtered signal
filtered_input_signal = sig.lfilter(FIR_coefficients, 1.0, input_signal)
# plot signals
if plot_time_signals :
x_ticks_range = np.arange(
0, (sample_nb+1)/sampling_rate, sample_nb/sampling_rate/10
)
y_ticks_range = np.arange(
2**(signal_bit_nb1), 2**(signal_bit_nb1)+1, 2**(signal_bit_nb4)
)
plt.figure("Time signals", figsize=(12, 9))
plt.subplot(2, 1, 1)
plt.title('input signal')
plt.step(time, input_signal, input_input_signal_color)
plt.xticks(x_ticks_range)
plt.yticks(y_ticks_range)
plt.grid()
plt.subplot(2, 1, 2)
plt.title('filtered signal')
plt.step(time, filtered_input_signal, filtered_input_signal_color)
plt.xticks(x_ticks_range)
plt.yticks(y_ticks_range)
plt.grid()
#
# Transfer function
#
# transfer function
(w, amplitude, phase) = transfer_function.bode(n=filter_order)
frequency = w / (2*np.pi)
# plot transfer function
if plot_transfer_function :
amplitude = np.clip(amplitude, 6*signal_bit_nb, 6)
log_range = np.arange(1, 10)
x_ticks_range = np.concatenate((1E2*log_range, 1E3*log_range, [2E4]))
amplitude_y_ticks_range = np.concatenate((
[6*signal_bit_nb],
np.arange(20*math.floor(6*signal_bit_nb/20), 1, 20)
))
phase_y_ticks_range = np.concatenate((
1E1*log_range, 1E2*log_range, 1E3*log_range
))
plt.figure("Transfer function", figsize=(12, 9))
plt.subplot(2, 1, 1)
plt.title('amplitude')
plt.semilogx(frequency, amplitude, transfer_function_amplitude_color)
plt.xticks(x_ticks_range)
plt.yticks(amplitude_y_ticks_range)
plt.grid()
plt.subplot(2, 1, 2)
plt.title('phase')
plt.loglog(frequency, phase, transfer_function_phase_color)
plt.xticks(x_ticks_range)
plt.yticks(phase_y_ticks_range)
plt.grid()
#
# Zeros
#
(zeroes, poles, gain) = sig.tf2zpk(FIR_coefficients, [1])
#print(zeroes)
# plot location of zeroes
if plot_zeros :
max_amplitude = 1.1 * max(abs(zeroes))
cutoff_angle = cutoff_frequency/sampling_rate*2*math.pi
cutoff_x = max_amplitude * math.cos(cutoff_angle)
cutoff_y = max_amplitude * math.sin(cutoff_angle)
plt.figure("Zeros", figsize=(9, 9))
plt.plot([max_amplitude, max_amplitude], [0, 0], locus_axes_color)
plt.plot([0, 0], [max_amplitude, max_amplitude], locus_axes_color)
plt.scatter(
np.real(zeroes), np.imag(zeroes),
marker='o', c=locus_zeroes_color)
plt.plot(
[cutoff_x, 0, cutoff_x], [cutoff_y, 0, cutoff_y],
locus_cutoff_frequency_color, linestyle='dashed'
)
plt.axis('square')
#
# Show plots
#
plt.show()
Filter structure edit
A Finite Impulse Response (FIR) filter's output is given by the following sequence:
The following figure shows the direct implementation of the above formula, for a 4^{th} order filter (N = 4):
With this structure, the filter coefficients are equal to the sample values of the impulse response.
IIR Filter Design
IIR filters are typically designed basing on continuoustime filter functions. Once the transfer function has been chosen, different filter structures allow to implement the filter, be it in hardware or in software.
Transfer function edit
The classical design technique is to start from a well known lowpass filter function and to transform it to a highpass, a bandpass or a bandstop function if required. Classical transfer functions are:
 Butterworth
 Chebyshev
 Elliptic or Cauer
The general form of the transfer function is given by the ratio of 2 polynoms:
Filter structures edit
The most straightforward digital filter implementation would be done by implementing its difference equation which is given by:
The following circuit shows a possible implementation of the difference equation for a 4^{th} order filter:
Nevertheless, this form is not used in practice. Indeed, the first who tried to implement a filter in that way have found out that this circuit is very sensitive to the coefficient values: a small change in one coefficient can have dramatic effects on the shape of the filter's transfer function. So rounding the coefficients to integer or fixedpoint values resulted into a design nightmare. Because of this, engineers turned over to one of the following approaches:
 simulating analog filter structures which had shown not to suffer of this acute sensitivity
 breaking the total transfer function into a chain of second order sections, with an additional first order section for evenorder filters
AllPole Filters edit
Allpole filters are limited to lowpass devices. The fact that they do not have poles makes their structure much simpler. Indeed, their transfer function is written as:
Some filter structures based on the simulation of analog filters are used for the implementation of allpole filters.
Butterworth and Chebyshev Type I functions are of allpole kind.
Chains of Integrators edit
A Chain of Integrators with Feedback (CIF) allows a straightforward implementation of an allpole transfer function.
This circuit, with to being the multiplier ouputs, can be described by the following equations:
which can be iteratively solved as:
From this, the time constants can be iteratively retrieved from the transfer function's coefficients.
The digital structure corresponding to the analog CIF is given in the following figure:
If the sampling rate is much higher than the time constants, the filter coefficients can be obtained directly from the CIF in the analog domain. Else, they can be found from the sampled system transfer function:
with .
Leapfrog Filters edit
A leapfrog filter simulates the functionality of an analog ladder filter and provides a robust implementation of an allpole transfer function. It can be shown that a large variation of coefficient values only leads to small variations of the transfer function, especially in the passband.
Taking the example of the following 4^{th} order allpole lowpass ladder filter,
one can write the following equations:
Rewriting these equations in order to have only voltages on the left side gives:
This new set of equations shows that the filter can be implemented with a chain of 4 integrators with:
 state variables having the form of or
 time constants having the form of or
The resulting circuit is an integrator chain where each integrator output is brought back to the preceeding integrator input. Hence the name leapfrog.
This analog device can be implemented as a digital circuit by replacing the analog integrators with accumulators:
The filter coefficients can be derived from tables providing the , and values for the selected filter functions and from the sampling period . An interesting point is that the relative values of the coefficients set the shape of the transfer function (Butterworth, Chebyshev, …), whereas their amplitudes set the cutoff frequency. Dividing all coefficients by a factor of two shifts the cutoff frequency down by one octave, which corresponds to a factor of two.
Let us note that the replacement of the analog integrators with accumulators corresponds is a very primitive design technique. In terms of signal processing, this corresponds to simplify the Ztransform to , which are the two first terms of the Taylor series of . Nevertheless, this approximation is good enough for filters where the sampling frequency is much higher than the signal bandwidth.
The state space representation of the preceeding filtre can be written as:
From this equation set, one can write the A, B, C, D matrices as:
From this representation, signal processing tools such as SciPy, Octave or Matlab allow to plot the filter's frequency response or to examine its zeroes and poles.
A special case is the Butterworth 3^{rd} order filter which has time constants with relative values of 1, 1/2 and 1. Due to that, this filter can be implemented in hardware without any multiplier, but using shifts instead.
Filters with poles and zeroes edit
Filters with poles and zeroes are best implemented in the form of a chain of biquadratic cells or by the parallel connection of two allpass filters.
Chain of Second Order Sections edit
A second order section, often referred as biquad, implements a second order transfer function. The transfer function of a filter can be split into a product of transfer functions each associated to a pair of poles and possibly a pair of zeroes. If the transfer function's order is odd, then a first order section has to be added to the chain. This section is associated to the real pole and to the real zero if there is one.
The second order section's transfer function is given by:
The most known biquad structures are shown next.
Directform 1 biquad edit
The following figure shows a direct form 1 biquad:
Directform 1 transposed biquad edit
The following figure shows a direct form 1 transposed biquad:
Directform 2 biquad edit
The following figure shows a direct form 2 biquad:
Directform 2 transposed biquad edit
The following figure shows a direct form 2 transposed biquad:
Obviously, the corresponding first order sections are done by removing the multipliers for coefficients and , together with the associated delay element(s).
Parallel connection of allpass filters edit
Allpass filters are used to modify the phase of a signal's frequency components without altering their amplitude. They can be used to compensate for other undesired phase shifts, as for example in the design of closetolinear phase IIR filters.
Additionally, the parallel connection of two allpass filters allows to implement any kind of transfer function. Indeed, at frequencies where the two branches show the same phase shift, the output will follow the input, but at the frequencies where the relative phase shift is 180°, the output will be zero. An example of this functionality is given in the family of Wave Digital Filters (WDFs).
Filter Representation
Direct Form I edit
Direct Form II edit
Transpose edit
Cascade Form edit
Polyphase Form edit
Parallel Form edit
Tunable Filters
We have highpass and lowpass filters, and it turns out that we can create both filter effects using a single piece of hardware called a tunable filter. Tunable filters rely on a device called an allpass filter to create other types of filter outputs.
Allpass Filters edit
A common, although unintuitive, tool in the DSP toolbox is the allpass filter. As the name implies, the allpass filter passes all frequency components without attenuation.
Why Allpass? edit
The question comes up, if an allpass filter doesn't actually filter anything, why do we bother to study it? The answer is that even though the allpass filter doesn't affect the magnitude response of the system, it will affect the phase response of the system.
Allpass filters have the property that cascading an allpass filter with another IIR filter will produce a linear output for the whole system.
Designing an Allpass edit
Tunable Filter Concept edit
Tunable Filter Design edit
Infinite Input Filtering
Filtering a data input that is infinite in length can be a tricky task for any computer to manage, and DSP engineers have identified two techniques for breaking down an infinite input into manageable chunks for filtering.
"Infinite Input" edit
infinite input response (IIR) filters are structures that use a feedback element to help filter data. An IIR filter has a transfer function of the form:
The constants a_{n} are known as the zeros of the transfer function, and the b_{m} terms are known as the poles of the transfer function. FIR filters, by comparison do not have poles.
IIR filters are named as such because the filter does not stop responding once the input has stopped, but instead the feedback element will continue to input values into the filter structure for processing after the input has stopped.
IIR vs FIR Filters edit
FIR  IIR  

Stability  Always BIBO Stable  May not be stable 
Efficiency  Not efficient  Very Efficient 
Phase  May have linear phase  phase not linear 
An FIR filter will have an order that is a significantly higher than the order of an equivalent IIR filter. This means that FIR filters need more complexity and more components then IIR filters to complete the same task. However, IIR filters have a potential for instability, which requires more analysis.
Overlap and Add edit
Overlap and Save edit
Comparison of Results edit
Ideal Filters
IDEAL DIGITAL SIGNALS edit
Filter Realizations edit
Realization Effects edit
Gibbs Phenomenon edit
Least Squares Design
The mathematical technique known as least squares can be used to design digital filters. This method works by creating the polynomial (transfer function) that is the "best fit" to a set of desired data (magnitude response).
Analog Filter Design
It may seem strange to have a section devoted to analog filter design in a digital signal processing book. However, there is a method to the madness. It turns out that digital filters can be modeled using analog filters, and analog techniques—which have been studied in depth for many years prior to digital systems—can be utilized to quickly and efficiently design digital filters.
The following chapters will show some of the techniques involved in designing analog filters, and transforming analog filters from one type to another. The section will culminate on a chapter on the bilinear transform, which transforms the mathematical representation for an analog filter into the equations for an equivalent digital filter.
Analog Design Process edit
Digital filters can be designed using analog design methods by following these steps:
 Filter specifications are specified in the digital domain. The filter type (highpass, lowpass,bandpass etc.) is specified.
 An equivalent lowpass filter is designed that meets these specifications.
 The analog lowpass filter is transformed using spectral transformations into the correct type of filter.
 The analog filter is transformed into a digital filter using a particular mapping.
There are many different types of spectral transformations and there are many mappings from analog to digital filters. the most famous mapping is known as the bilinear transform, and we will discuss that in a different chapter.
Butterworth edit
Butterworth ensures a flat response in the passband and an adequate rate of rolloff. A good "all rounder," the Butterworth filter is simple to understand and suitable for applications such as audio processing.
Chebyshev I edit
The Chebyshev 1 filter has ripple in the passband of the filter.
Chebyshev II edit
Chebyshev 2 has ripple in the stopband.
Elliptic edit
This filter has equiripple (the same amount of ripple in the passband and stopband).
Further reading edit
Spectral Transforms
When designing a particular filter, it is common practice to first design a Lowpass filter, and then using a spectral transform to convert that lowpass filter equation into the equation of a different type of filter. This is done because many common values for butterworth, cheybyshev and elliptical lowpass filters are already extensively tabulated, and designing a filter in this manner rarely requires more effort then looking values up in a table, and applying the correct transformations. This page will enumerate some of the common transformations for filter design.
It is important to note that spectral transformations are not exclusively used for analog filters. There are digital variants of these transformations that can be applied directly to digital filters to transform them into a different type. This page will only talk about the analog filter transforms, and the digital filter transforms will be discussed elsewhere.
LowPass to LowPass edit
This spectral transform is used to convert between two lowpass filters with different cutoff frequencies.
LowPass to HighPass edit
LowPass to BandPass edit
LowPass to BandStop edit
Analog Mappings
There are a number of different mappings from the analog domain to the digital domain, and vice versa. All mappings must satisfy the following rules:
 Stability must be preserved in both domains. This means that the interior of the unit circle in the digital domain must correspond to the lefthalf plane in the analog plane.
 The mapping must be onetoone. every point in one domain must correspond to a unique point in the other domain. This relationship should hold both ways.
The relationship between the analog and digital domains is highly nonlinear, and there is no way to perfectly reproduce an analog filter in a digital system, or to reproduce a digital filter in an analog system. However, it is possible to faithfully reproduce (within a certain tolerance) the magnitude response of the filter. It is rarely possible to reproduce the phase response of a filter across the transformation, however.
Bilinear Transform
The Bilinear Transform edit
The Bilinear transform is a mathematical relationship which can be used to convert the transfer function of a particular filter in the complex Laplace domain into the zdomain, and viceversa. The resulting filter will have the same characteristics of the original filter, but can be implemented using different techniques. The Laplace Domain is better suited for designing analog filter components, while the ZTransform is better suited for designing digital filter components.
The bilinear transform is the result of a numerical integration of the analog transfer function into the digital domain. We can define the bilinear transform as:
The bilinear transform can be used to produce a piecewise constant magnitude response that approximates the magnitude response of an equivalent analog filter. The bilinear transform does not faithfully reproduce the analog filters phase response, however.
As example, the phase response for the bilinear equivalent to an analog system is shown next, using a sampling frequency of 40 rad/sec. Compare the phase responses in the plot, realizing the limit is Ws/2=20 rad/s  the clear divergence is evident.
 num = [1 11.21 116.242 372.601 561.589 363.528 ];
 den= [1 26.489 340.47 2461.61 10433.1 23363.9 19049. 4981.82 ];
 C=tf(num,den); % GNU Octave command
 D=c2d(C,2*pi/40,'bi'); % A/D via Octave
Note that two zeros have been added so that zero and pole counts match.

Phase discrep approaching Ws/2
W Domain edit
The bilinear transform is a onetoone mapping, that is that a unique point in one domain will be transformed into a unique point in the other domain. However, the transformation is not a linear transformation, and is not an exact equivalency between Laplace and Z domains. If a digital filter is designed and transformed into an analog filter, the resulting poles and zeros are not located on the splane, but instead on the wplane that has similar properties, but with a nonlinear correspondence.
In the bilinear transform, the positive imaginary axis in the sdomain is transformed into the upperhalf unit circle in the zdomain. Likewise, the negative imaginary axis in the sdomain is transformed into the lowerhalf unit circle in the zdomain. This mapping is highly nonlinear, however, resulting in a phenomena known as "frequency warping".
Prewarping edit
Frequency warping follows a known pattern, and there is a known relationship between the warped frequency and the known frequency. We can use a technique called prewarping to account for the nonlinearity, and produce a more faithful mapping.
The p subscript denotes the prewarped version of the same frequency. Note that in the limit , the continuous solution is .
Phase edit
The bilinear transform does not maintain phase characteristics of the analog filter, and there is no way to correct the phase response to match.
Filter Design Steps edit
When designing a digital filter using an analog approximation and the bilinear transform, we follow these steps:
 Prewarp the cutoff frequencies
 Design the necessary analog filter
 Apply the bilinear transform to the transfer function
 Normalize the resultant transfer function to be monotonic and have a unity passband gain (0dB).
Alternatively, if we have an inverse bilinear transform, we can follow these steps:
 Use the inverse bilinear transform on the filter specifications in the digital domain to produce equivalent specifications in the analog domain.
 Construct the analog filter transfer functions to meet those specifications.
 Use the bilinear transform to convert the resultant analog filter into a digital filter.
The Inverse Bilinear Transform edit
The inverse bilinear transform can be specified as such:
Where is the sample rate
This can also be written as:
which very much resembles the Bilinear Transform.
Discrete Cosine Transform
Discrete Cosine Transform in dsp edit
The Discrete Cosine Transform (DCT) is a transform that is very common when encoding video and audio tracks on computers. Many "codecs" for movies rely on DCT concepts for compressing and encoding video files. The DCT can also be used to analyze the spectral components of images as well. The DCT is very similar to the DFT, except the output values are all real numbers, and the output vector is approximately twice as long as the DFT output. It expresses a sequence of finite data points in terms of sum of cosine functions.
Inverse DCT edit
Computing DCT edit
Uses of DCT edit
The JPEG process is a widely used form of lossy image compression that centers around the Discrete Cosine Transform. The DCT works by separating images into parts of differing frequencies.
Further reading edit
Haar Transform
Haar Wavelet edit
Haar Matrix edit
Haar Transform edit
The Haar Transform, or the Haar Wavelet Transform (HWT) is one of a group of related transforms known as the Discrete Wavelet Transforms (DWT). DWT Transforms, and the Haar transform in particular can frequently be made very fast using matrix calculations. The fastest known algorithm for computing the HWT is known as the Fast Haar Transform, and is comparable in speed and properties to the Fast Fourier Transform.
Uses of the Haar Transform edit
Haar transform uses nonsinusoidal basic wavefunction. So it has great applications related to DSP. The basic haar transform matrix is defined by the function Hk(x). Where o<=k<=N1, N is the matrix size.
Computing Haar Transform edit
Further reading edit
Quantization
FixedPoint Numbers edit
Floating Point Numbers edit
Dynamic Range Scaling edit
Quantization Errors edit
Further reading edit
 Wikibooks: Floating Point
 Wikibooks: Embedded Systems/Floating Point Unit
Multirate DSP
Rate Alteration edit
Arbitrary Rate Circuits edit
Nyquist Filters edit
QuadratureMirror Filters edit
LChannel Filters edit
Windowing
Most digital signals are infinite, or sufficiently large that the dataset cannot be manipulated as a whole. Sufficiently large signals are also difficult to analyze statistically, because statistical calculations require all points to be available for analysis. In order to avoid these problems, engineers typically analyze small subsets of the total data, through a process called windowing.
Windowing Introduction edit
Windowing is the process of taking a small subset of a larger dataset, for processing and analysis. A naive approach, the rectangular window, involves simply truncating the dataset before and after the window, while not modifying the contents of the window at all. However, as we will see, this is a poor method of windowing and causes power leakage.
Applying Windows edit
Application of a window to a dataset will alter the spectral properties of that dataset. In a rectangular window, for instance, all the data points outside the window are truncated and therefore assumed to be zero. The cutoff points at the ends of the sample will introduce highfrequency components.
Consider the system H(z), with input X(z) and output Y(z). We model this as:
If we have a window with transfer function W(z), we can mathematically apply the window to our signal, X(z) as such:
Then, we can pass our windowed signal into our system, H(z) as usual:
Leakage edit
If our signal is a lowpass or passband signal, the application of a window will introduce highfrequency components. Power from the original signal will be diverted from the specified frequency band into the highfrequency areas. This redistribution of power from a target band to the upper frequencies is known as leakage.
If we look at a rectangular window, we know from duality that the frequency response of that window is a sinc function, which has nonzero values all the way out to both positive and negative infinity. Convolution of the sinc function with any narrowband signal is going to cause a very spreadout spectrum.
Types of Windows edit
Hamming edit
function of hamming
0.54+0.46cos2pi n/m1
Blackman edit
0.42+.5cos(2*pi*n/m1)+.08cos(4*pi*n/m1)
Greatest stop band attenuation of mentioned windowing techniques at the expense of a larger transition band.
The eq. for the symmetric 4term Blackman Harris window of length N is....
w(n) = a0  a1*cos(2*pi*n/N1) + a2*cos(4*pi*n/N1)  a3*cos(2*pi*n/N1); 0 <= n <= N1
a0 = 0.35875 a1 = 0.48829 a2 = 0.14128 a3 = 0.01168
SigmaDelta modulation
A SigmaDelta Modulator (ΣΔ modulator) allows to operate Analog to Digital Conversion (ADC) or Digital to Analog Conversion (DAC) by the means of a onebit signal.
The usage of a single bit signal is also used by PulseWidth Modulation (PWM), where the SignaltoNoise Ratio (SNR) is worse but the switching slower.
First Order Modulator edit
Circuit edit
The first order ΣΔ modulator consists out of an accumulator and a comparator.
The output is a bit stream and the original signal can be reconstruted by the means of a lowpass filter.
Analysis edit
The frequency response of the sigmadelta modulator can be analysed by replacing the comparator with the addition of an error signal.
From these equations, one finds the Signal Transfer Funtion (STF)
and the Noise Transfer Funtion (NTF)
The STF is a one period delay but the NTF is a first order highpass function. This indicates that the modulation noise, due to the switching, occupies higher frequencies. If the input signal is limited to low frequencies, it is possible to separate the signal from the nois using a lowpass filter. As the highpass transfer function decreases by 20 dB/decade (or 6 dB/octave) as the frequency decreases, the signal to noise ratio increases by the same rate. Integrating the noise amplitude over the signal band shows that a first order modulator gains a signal to noise ration of 1.5 bit per octave. This can be compared to PCM which shows a gain of only 1 bit per octave.
When its input is constant, a first order modulator will provide a cyclic pattern. This pattern can become quite long for specific input values, and this leads to lowfrequency ringing tones in the modulated signal. As these tones go lower in frequency, they become more and more difficult to separate from the original signal. Higher order modulators show less repetitive patterns and are thus preferred to first order ones.
Second Order Modulator edit
Circuit edit
A second order ΣΔ modulator will require 2 accumulators. Different topologies are possible. The following circuit shows a typical second order circuit.
The 2 coefficients allow to control the digital to analog transfer function together with the corresponding noise shaping characteristic.
Analysis edit
Here too, the comparator can be replaced by the insertion of an error signal:
The NTF and STF can be found from the following equations:
which can be rewritten as:
Signal Transfer Function edit
The STF is found by writing , which gives:
This can be rewritten in the state space representation as:
The statespace matrices are built from the modulator matrices as:
Solving the system for the second order modulator gives:
Noise Transfer Function edit
The NTF is found by writing , which gives:
This can be rewritten in the state space representation as:
The statespace matrices are built from the modulator matrices as:
Solving the system for the second order modulator gives:
Comments edit
The matrices contain all the information to both simulate and analyse the behaviour of any sigmadelta modulator. The derivation of the statespace description of the NTF and the STF from these matrices is generally applicable to any sigmadelta modulator.
The matrices and are identical, which means that the STF and the NTF share the same poles. The STF has only poles and this results in a lowpass function. The NTF has 2 poles at , which gives a high pass response. This ensures that the signal and the noise are well separable.
The STF has a DC gain of:
A STF DC gain of 1 is achieved by multiplying the input by .
Design edit
The poles of the STF are given by
This equation allows to arbitrarily place the poles of the STF. For a pair of poles at:
the modulator coefficients are:
As an example, an allpole Butterworth STF with cutoff frequency at 1/4 the sampling frequency give poles:
and coefficients:
An allpole close to Bessel STF with cutoff frequency around 1/4 the sampling frequency give coefficients:
Note that the usual habits of using coefficients or both result in having poles exactly on the unit circle, which is not recommended in terms of stability.
Higher order topologies edit
For higher order modulators, different topologies exist^{[1]}.
Simplified Chain of Integrators with FeedBack (SCIFB) edit
The following picture shows a Simplified Chain of Integrators with FeedBack (SCIFB) modulator ^{[2]}:
With coefficient different from zero, this is a structure with a local resonator.
For this structure, the system equations are:
and the matrices are:
References edit
 ↑ Norsworthy, Steven R.; Schreier, Richard; Temes, Gabor C. (1997). DeltaSigma Data Converters. IEEE Press. ISBN 0780310454.
 ↑ Liu, Mingliang (2003), The Design of DeltaSigma Modulators for MultiStandard RF Receivers
Multirate Filters
Multirate Filters edit
There are many instances where the rate at which a signal can be processed by a particular component or module is different from the speed at which data enters that module.
Example: DVD and CD Audio
CD audio is encoded at 44.1K samples per second. DVD audio, by comparison, is encoded at 48K samples per second. In order to convert music from one to the other, the sample rate needs to be altered.
Upsampling edit
Upsampling creates more samples in the same amount of time, typically by inserting zerovalued samples between the preexisting samples. At first blush, this may seem like an illconceived approach. However, if one considers that a discrete signal is already zero between the sample points, the approach begins to make more sense.
Spectral Effects edit
Inserting a single zero between each of the samples will cause the spectrum to replicate and fold, creating a mirror image. If the original sample rate was , the new sample rate will be , and the spectrum will fold and replicate at . This image can then be filtered out using a lowpass filter.
Inserting more than one zero between sample points will cause the spectrum to replicate and fold accordionstyle. If N zeros are inserted, N images will be created. Again, these images can be attenuated by applying a lowpass filter.
Implementation edit
Typically, the lowpass filter and the upsampler are implemented as a unit, and the upsampling takes place only conceptually. Since the zerovalued samples will not contribute to the sum of products in an FIR filter, these multiplications are simply skipped. Zero times anything is zero.
Downsampling edit
Downsampling, or decimation is the process of discarding certain samples so that there are fewer samples in the same amount of time. Once discarded, those samples can never be replaced and error is introduced into the system. However, a downsampled system can also be processed with a slower filter, which is typically less expensive.
Spectral Effects edit
Downsampling causes the spectrum to spread. If the spectrum is periodic, there could be some overlapping of spectral objects, and this causes aliasing. Aliasing of this sort is typically resolved by passing the downsampled signal through a lowpass filter to help remove the overlapped areas.
Ideal Reconstructor edit
An ideal reconstructor can be created by having an upsampler followed directly by a downsampler.
Similarly, if a single is separated, each branch can be delayed and downsampled, and then combined together with zero loss. This enables the signal to be processed at a much slower rate then the input data rate, without losing any data.
Fractional Decimation edit
Upsampling and downsamping alter the size of the data set by an integer ratio of samples. In order to achieve a fractional sample rate, upsamplers and downsamplers need to be coupled together to change the data rate to a fraction of the input data rate.
Wiener Filters
A Wiener Filter is a filtering system that is an optimal solution to the statistical filtering problem.
Statistical Filtering edit
Statement of Problem edit
The operator E[] is the expectation operator, and is defined as:
where f_{x}[n] is the probability distribution function of x.
d[n] is the expected response value, or the value that we would like the input to approach.
e[n] is the estimation error, or the difference between the expected signal d[n] and the output of the FIR filter. We denote the FIR filter output with a hat:
Where the convolution operation applies the input signal, u[n], to the filter with impulse response w[n].
We can define a performance index J[w] which is a function of the tap weights of the FIR filter, w[n], and can be used to show how close the filter is to reaching the desired output. We define the performance index as:
J[w] is also known as the meansquared error signal. The goal of a Wiener filter is to minimize J[w] so that the filter operates with the least error.
Adaptive Filtering edit
Adaptive filtering is a situation where the coefficients of a filter change over time, typically in response to changes in the characteristics of the input signal.
WienerHopf Equations edit
Where R is the autocorrelation matrix. w_{o} is the optimal set of filter coefficients, arranged in a vector, and p is the crosscorrelation vector.
Wiener Filters edit
Wiener Filters are typically implemented with FIR filter constructions. This is because the wiener filter coefficients change over time, and IIR filter can become unstable for certain coefficient values. To prevent this instability, we typically construct adaptive filters with FIR structures.
DSP Programming
Digital Signal Processing units, on first glance, are very similar to other programmable microprocessors, microcontrollers, and FPGAs. However, DSP chips frequently have certain features and limitations involved that other categories of chips don't have. This chapter will attempt to explain, in broad languageneutral terms, some of the issues involved with DSP programming.
Saturation Arithmetic edit
Intelcompatible microprocessors, something that most programmers are familiar with, have several features that people have come to expect from programmable chips. In particular, Intelcompatible chips have an arithmetic mode known as "RollOver Arithmetic", also called modular arithmetic. Let's say that we have the integer number 250 (decimal) stored in a bytesized register. If we add 10 to that number, the processor will overflow, and will rollover the answer, giving us the answer 4. However, DSP chips frequently do not rollover, and instead saturate to the maximum value.
Let's say we are processing a grayscale image. This image has bytesized integer values of white/black value as such: Pure white is 255, and pure black is 0. Now, let's say we want to "brighten" every pixel by 10 points, to improve visibility. So we create a loop, and we go through each pixel in the image, adding 10 to the value at any particular point. But what happens if we add 10 to a pixel that is already pure white? In a regular Intelcompatible chip, the value 255 + 10 = 9. In effect, if we make the white pixels brighter, we are turning them black! now, to avoid this, DSP chips saturate, and will go up to the highest value (255), and won't roll over. So, on a DSP chip, white can't get any whiter. Here are some examples (using unsigned bytes):
Multiply and Accumulate (MAC) edit
In difference equations, we have seen that multiply and addition operations are the most common. Let's look at a general example:
Now, we notice that each element in this equation is multiplied by a coefficient, and added to the total. DSP engineers have noticed this pattern, and have optimized DSP chips to be able multiply and add very quickly (at the expense of other operations missing or slow). In fact, many DSP chips will have instructions known as multiply and accumulate, which will perform both operations simultaneously, much quicker than a normal processor could do.
Sound Processing
Digital Sound edit
Sampling Frequency edit
Sound in the digital realm is stored in one or more arrays of discrete samples, with each array of samples correlating to a channel (e.g. stereo sound requires two channels, and thus two arrays of samples). The interval of time between each sample is a constant, and is determined by the type of data to be represented. Since we are interested in sound, and the extreme upper limit of human hearing is generally accepted as 20 kHz, the NyquistShannon sampling theorem can be used to determine the interval between samples to accurately reconstruct the signals we're interested in.
This theorem states that
“  Exact reconstruction of a continuoustime baseband signal from its samples is possible if the signal is bandlimited and the sampling frequency is greater than twice the signal bandwidth.  ” 
Essentially what this means is a signal that is limited to a certain range (audible sound: ~20 Hz to 20 kHz) can be reconstructed without error if it is sampled at a rate that is greater than twice the bandwidth. The Red Book audio CD standard sets the sampling rate at 44,100 Hz. This frequency was chosen to leave ample overhead (as required by the NyquistShannon theorem), but could support at least up to 22 kHz.
44.1 kHz is the general standard for sampling rates in digital audio on consumer level equipment, however 48 kHz is common when working with film or video. Also, many recording engineers prefer to record classical or otherwise complex music at 88.2 or 96 kHz—some claim to be able to perceive a difference.
When converting from 48 kHz to 44.1 kHz a sonic blurring effect can sometime occur, because the math is floating point, which is inherently imprecise on a computer. The conversion from 88.2 kHz to 44.1 kHz or 96 kHz to 48 kHz is simpler to perform, since the computer, or device, doing the conversion only has to disregard half the samples. To bypass this problem, a highquality digitaltoanalog converter can be used to bring, for example, a 48 kHz signal back to its analog form, and then is fed into another highquality analogtodigital converter to resample the signal at 44.1 kHz. This technique is common practice in recording studios where highend equipment can be trusted to do the conversion flawlessly, however in other situations, the sonic distortion caused by converting the audio in software or hardware may be of little concern.
Bits Per Sample edit
While sampling frequency determines the time component of an audio signal, the number of bits per sample is used to describe the amplitude. Red Book audio CDs store each sample as a 16 bit signed integer. This means that when an audio signal is converted for use on a CD, each sample's value is quantized as an integer to fit in the range 32768 to +32767.
Wave Files edit
Wave files contain data which is a representation of audio sound. This format for storing data is an uncompressed format. This means the data can be sent to the digitaltoanalog processor for playback without an added step of decompression. This also means that this format will consume a great deal of memory.
MP3 Compression edit
OGG Compression edit
Image Processing
Image Processing Techniques edit
BMP Images edit
JPEG Images edit
GIF Images edit
PNG Images edit
Digital Signal Processors
Special DSP Hardware edit
Saturation Arithmetic edit
Multiply and Accumulate edit
Specific Models edit
Transforms
This page lists some of the transforms from the book, explains their uses, and lists some transform pairs of common functions.
ContinuousTime Fourier Transform (CTFT) edit
[CTFT]
CTFT Table edit
Time Domain  Frequency Domain  

1  
2  
3  
4  
5  
6  
7  
8  
9  
10  
11  
12  
13  
14  
15  
16  
Notes: 
 
DiscreteTime Fourier Transform (DTFT) edit
DTFT Table edit
Time domain where 
Frequency domain where 
Remarks 

Definition  
Here represents the delta function which is 1 if and zero otherwise.  
is 2π periodic  
DTFT Properties edit
Property  Time domain 
Frequency domain 
Remarks 

Linearity  
Shift in time  integer k  
Shift in frequency  real number a  
Time reversal  
Time conjugation  
Time reversal & conjugation  
Derivative in frequency  
Integral in frequency  
Convolve in time  
Multiply in time 