Lab Man DSP

Download as pdf or txt
Download as pdf or txt
You are on page 1of 89

Digital Signal Processing Lab Manual

Contents
List of figure:- ........................................................................................................................................ 2
List of Table:- ........................................................................................................................................ 4
Tutorial 1: Introduction to MATLAB ........................................................................................................ 5
Tutorial 2: Introduction to MATLAB ................................................................................................ 13
Experiment Part one: Generation of Discrete-Time Signals in Time Domain ........................................ 26
Lab #1: Generation of Discrete-Time Signals in Time Domain ........................................................... 26
Experiment Part Two: Simulation of Discrete-Time Systems............................................................ 36
Lab#1The Moving Average System ................................................................................................ 36
Lab #2 Discrete-Time Systems in the Time Domain for LTI ................................................................ 39
Lab #3 Linear Convolution ............................................................................................................. 44
Lab #3:- Circular Convolution ........................................................................................................ 46
Experiment Part Three: The Frequency Domain Analysis ................................................................ 49
Lab #1: Discrete-Time Fourier Transform (DTFT) Computation ......................................................... 49
Lab #2:-The Discrete Fourier Transform (DFT) and Inverse DFT Computations ................................. 53
Lab #4:- Z-Transform Analysis ...................................................................................................... 57
Experiment Part Four: Sampling Process .......................................................................................... 61
Lab #1Sampling of Continuous-Time Signal into Discrete-Time Signal .............................................. 63
Lab #2 Aliasing Effect in Time Domain .......................................................................................... 65
Lab #3 Effect of Sampling in Frequency Domain: Aliasing Effect in Frequency .......................... 67
Lab #4 DECIMATION & INTERPOLATION .............................................................................. 70
Experiment Part Five: Digital Filter Design ............................................................................................ 76
Lab #1 Infinite Impulse Response (IIR) Filter Design Based on the Impulse Invariance Method
and Bilinear Transformation .......................................................................................................... 76
Lab #2 Filter Design by Matlab tools ................................................................................................. 84

1 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

List of figure:-
Figure T1. 1 Window of the MATLAB 6
Figure T1. 2 graphical help window of the matlab 7
Figure T2. 1 Graphical presentation using plot(x,y) function of the MATLAB 19

Figure T2. 2 Graphical presentation using plot(x,y,’r*’) function of the MATLAB 19

Figure T2. 3 Graphical presentation using stem(x,y) function of the MATLAB. 19

Figure T2. 4 Graphical presentation by plot function of the MATLAB. 20

Figure T2. 5 Graphical presentation by plot function of the MATLAB. 21

Figure T2. 6 Graphical presentation by plot function of the MATLAB. 22

Figure T2. 7 Edit Window with a program to demonstrate the animated plot 23

Figure T2. 8 Animated plot of m-file tutoria1212.m 23

Figure T2. 9 Animated plot using command ‘comet3’ 24

Figure P1L1. 1 Impulse response-sample output of program1.1 28

Figure P1L1. 2 Impulse response-sample output of program1.2 29

Figure P1L1. 3 The unit step response: sample output of program 1.3 30

Figure P1L1. 4 Discrete-time representation of the unit ramp signal. 31

Figure P1L1. 5 Complex exponential sequence: sample plot of program 1.5. 32

Figure P1L1. 6 Sinusoidal sequence: sample plot of program 1.6. 33

Figure P1L1. 7 Sinusoidal sequence: sample plot of program 1.7 34

Figure P1L1. 8 square wave 35

Figure P1L1. 9 Sawtooth Wave Sequence 36

Figure P2L1. 1 The illustration of the moving average filter Program 2.1 . 38

Figure P2L1. 2 Demonstration of linearity property of the discrete-time system . 41


Figure P2L1. 3 Demonstration of Time-Invariant and Time-Varying Systems of the discrete-time system
42

Figure P2L1. 4 Demonstration Stability of LTI Systems of the discrete-time system. 43

Figure P2L1. 5 linearly convolved output of LTI system to a finite length input. 45

Figure P2L1. 6 Convoluted sum by using the circular convolution. 48

Figure P3: 1 Compute the frequency samples of the DTFT 51

2 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Figure P3: 2 The DFT of the input sequence x(n) 54

Figure P3: 3 FFT using DSP function fft 55

Figure P3: 4 compute FFT 56

Figure P3: 5 Pole- zero diagram of H (z) 59

Figure P4: 1 Sampling of continuous time signal to discrete-time signal. 64

Figure P4: 2 Aliasing Effect in the Time Domain 65

Figure P4: 3 Effect of aliasing in frequency domain 68

Figure P4: 4 Illustration of the Sampling Process 71

Figure P4: 5 Illustration of Down-Sampling by an Integer Factor 71

Figure P4: 6 Illustration of Decimation Process 72

Figure P4: 7 Illustration of Interpolation Process 73

Figure P4: 8 Illustration of Sampling Rate Alteration 74


Figure P5. 1 Magnitude responses of IIR filter using impulse invariance and bilinear transformation by
method I. 81
Figure P5. 2 The magnitude responses of IIR filter using impulse invariance and bilinear transformation
by direct method 82

Figure P5. 3 SPTool startup 85

Figure P5. 4 Filter Designer 85

Figure P5. 5 Designing filter 86

Figure P5. 6 show low pass filter 87

Figure P5. 7 Export way87

Figure P5. 8 show a low pass filter 88

3 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

List of Table:-
Table P5: 1 Analogue Transformation Functions ................................................................................... 78
Table P5: 2 Transformation from Analogue to Digital Filters ................................................................. 79
Table P5: 3 Direct Design Matlab Functions for IIR Filters .................................................................... 79

4 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Tutorial 1: Introduction to MATLAB


Objective & Aim
 To familiarize and practicing of the general MATLAB functions.

T1.1 Introduction to MATLAB


MATLAB is a package that has been purpose-designed to make computations easy, fast and
reliable. The idea behind MATLAB was to provide a simple way of using these programs that
hide many of the complications. The idea was appealing to scientists and engineers who needed
to use high performance software but had neither the time nor the inclination (nor in some cases
the ability) to write it from scratch. Since its introduction, MATLAB has expanded to cover a
very wide range of applications. Now, it can be used as a very simple and transparent
programming language where each line of code looks similar to the mathematical statement it is
designed to implement.
Basic MATLAB is appropriate for the following.
 Computations, including linear algebra, data analysis, signal processing, polynomials and
interpolation, numerical integration (quadrature), and numerical solution of differential
equations.
 Graphics, in 2-D and 3-D, including color, lighting, and animation.
It also has collections of specialized functions, called toolboxes that extend its functionality. In
particular, it can do symbolic algebra, e.g. it can tell you that (x+y)^2 is equal to
x^2+2*x*y+y^2.
The purpose of this course is to make to understand the DSP deeply. Therefore, we won’t be
using the built-in functions unless it is necessary. MATLAB programs for complicated tasks may
be somewhat slower to run than programs written in languages such as C. However, the
MATLAB programs are very easy to write. We shall use MATLAB as a vehicle for learning
elementary programming skills and applications.
These are skills which will be useful independently of the language you choose to work in. Some
students may already know some of these skills, but we shall not assume any prior knowledge.
MATLAB has four main windows.
1. Command window.
2. Edit or Debug window.
3. Graphical or Figure window.
4.Current Folder or Work Space

5 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

However, the Command Window has several other windows, which is shown in Figure T1. 1.
Both Edit Window and Graphical Window modify the existing ones. Similarly, a
Figure/Graphical Window is used to display MATLAB graphics.
The Command Window is the main window of MATLAB. A user can enter MATLAB
commands at the command prompt (>>) in the Command Window and they will be executed. An
Edit/Debug Window is used to create a new MATLAB files (m-files) or to

Editor
Current
Folder

Command
Window
Work
spaceFigure T1. 1 Window of the MATLAB.

To go through MATLAB, one doesn’t need a book. The help facilities available inside the
MATLAB itself is more than sufficient. MATLAB has two ways of getting help; one graphical
help and the next is text help. The graphical help is available from the Command Window itself.
However, the text help is obtained from the command prompt (>>). The Graphical Help Window
is shown in Figure T1. 3.

6 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Figure T1. 2 graphical help window of the matlab.


Some of the most useful MATLAB commands are given below:
On-line help
 Help……………… lists topics on which help is available
 Helpwin……………..opens the interactive help window
 Helpdesk…………… opens the web-browser-based help facility
 Help…………………topic provides help on topic look for string lists
Workspace information and control
 Who…………………… lists variables currently in the workspace
 Whos…………………… as above, giving their size
 What……………………. lists m-files on the disk
 Clear……………………. clears the workspace, removing all variables
 clear all…………………. clears all variables and functions from the workspace
Command window control
 clc………………………….. clears command window, command history is lost
 home……………………….. same as clc
Graphics
 plot…………………………….. plots a graph
 stem…………………………… plots a discrete graph
 xlabel('x')……………………… labels x axis x
 ylabel('y')……………………… labels y axis y
 title('title')……………………… gives a figure a title
 axis………………………. fixes figure axes
 clf………………………… clears figure from figure window

7 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

 cla………………………… clears figure from figure window, leaving axes


Controlling program execution
 break…………………………………… terminates execution of a for or while loop
 error('message')………………………… aborts execution, displays message on screen
 return…………………………………… exit from function, return to invoking program
Input from and output to terminal
 x=input('type x:')…………………………….. asks user to give a value to be assigned to
x
 disp('string')…………………………………. outputs string to terminal
String-number conversion
 num2str……………………………………… converts a number to a string (so it can be
output e.g. as part of a message)Logical functions
 isempty true (=1)………………………….. if a matrix is empty find finds indices of non-
zero elements of a matrix
Arithmetic functions
 sum(x)……………………………….. calculates the sum of the elements of the vector x
 prod(x)……………………………... calculates the product of the elements of the vector
x
Termination
 ^c (Control-c)……………………….. local abort, kills the current command execution
 Quit………………………………….. quits MATLAB
 Exit…………………………………… same as quit

T1.2 Familiarization with MATLAB Functions


 To familiarize with MATLAB, let us go through Tutorials.

T1.2.1 Tutorial 1: Simple Calculations


In this tutorial, you will start to become familiar with the MATLAB development environment
and some of its facilities. You will learn
 how to start and quit MATLAB,
 how to do simple arithmetic calculations,
 how to perform common mathematical functions,
 how to assign values to variables, and
 how to use some of MATLAB's built-in functions.

8 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

T1.2.1.1 Starting and Closing MATLAB


From start menu, find the MATLAB icon and click it. The MATLAB Command Window
with a prompt (>>), will be displayed. That is, your MATLAB is ready to use. Please confirm
that only one Command Window has been opened. Do not open two or more than two Command
Windows. To finish MATLAB session, type exit in the prompt (>>) and press enter. Instead, one
can close the window directly. However, exit is the safest way of closing the MATLAB window.
From now on, when you are instructed to type various commands, unless you are explicitly
told to type them elsewhere, they should be typed to the prompt >> in the MATLAB Command
Window.
T1.2.1.2 Simple Arithmetic in MATLAB
The basic operations are +, -, *, /, and ^, which stand for add, subtract, multiply, divide and
exponentiate, or raise to the power of. Work through the following examples, which show the
results of typing some simple arithmetic commands to the MATLAB prompt. The commands
typed by the user are those immediately following the prompt, >>. A % symbol means a
comment: the rest of this line is ignored by MATLAB. When a computation produces a
response, this is displayed immediately below the command that produced it. Here are some
examples:
>> clear % removes any old variables from the workspace
>> format short % outputs results in decimal form with 5 digit accuracy
>> 2+2 % two numbers get added and result displayed on screen
ans =
4
>> 4-5 % subtracted and result displayed on screen
ans =
-1
>> 4*5 % multiplied and result displayed on screen
ans =
20
>> 4/5 % divided and result displayed on screen
ans =
0.8000
>> 3^2 % exponentiated and result displayed on screen
ans =
9
>> (2+3i)^2 % exponential to complex and result displayed on screen
ans =
-5.0000 +12.0000i
FORMAT Set output format
All computations in MATLAB are done in double precision. FORMAT may be used to switch
between different output displays formats as follows:
 format Default. Same as SHORT.

9 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

 format short Scaled fixed point format with 5 digits.


 format long Scaled fixed point format with 15 digits.
 format short e Floating point format with 5 digits.
 format long e Floating point format with 15 digits.
 format short g Best of fixed or floating point format with 5 digits.
 format long g Best of fixed or floating point format with 15 digits.
 format hex Hexadecimal format.
 format bank Fixed format for dollars and cents.
 format rat Approximation by ratio of small integers.
Spacing:
 format compact Suppress extra line-feeds.
 Format loose Puts the extra line-feeds back in.

T1.2.1.3 Common Mathematical Functions in MATLAB


MATLAB contains a large number of mathematical functions. Most of them are entered exactly
as you would write them mathematically. For example,
>> sin(3) % sine function, result displayed in screen
ans =
0.1411
>> exp(2) % exponential function, result displayed in screen
ans =
7.3891
>> log(10) % logarithmic function, result displayed in
screen
ans =
2.3026

T1.2.1.4 Assign Values to Variable in MATLAB


Variables can be used to store numerical values and matrices in MATLAB. Here are some
examples:
>> x = 2+2 % result loaded into new variable x and displayed
x=
4
>> y = 2^2 + log(pi)*sin(x); % more complicated computation illustrates

10 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

% built-in functions log and sin and built-in


% constant pi. The; at the end suppresses output
>> y % prints the contents of y
y = 3.1337
>> format short e % changes format of output
>> y
y = 3.1337e+00
>> format long % changes format of output again
>> y
y = 3.13366556593561
>>more on % this allows text to be presented page by page on the
% screen. Without it the text just scrolls around very
% quickly and is impossible to read
>> help format % prints the contents of the manual for the
% command ``format''
>> who % displays the current variables in the workspace
>> whos
Name Size Bytes Class Attributes

ans 1x1 16 double complex

T1.2.1.5 Some of MATLAB’s Build-in Functions


Type the following functions in the Command Window and see the responses:
>>help
>>help plot
>>help colon
>>help ops

11 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

>>help punct
>>help zeros
>>help ones

T1.2 Exercises for Tutorial 1


1. Using MATLAB, do the following calculations:
(i) 1 + 2 + 3; (ii) 2 ; (iii) cos(π 6)
Note that MATLAB has pi as a built-in approximation for π (accurate to about 13 decimal
places), and built-in functions sqrt and cos. Look up the manual pages for these if necessary.
2. After MATLAB has done a calculation the result is stored in the variable ans. After you have
done (iii) above, type ans^2 to find cos2 (π 6) . The true value of this is 3/4 . (Check that you
know why!) Print out ans using both format short and format long and compare with the true
value. On the other hand use MATLAB to compute ans-3/4. This illustrates rounding error: the
fact that the machine is finite means that calculations not involving integers are done inexactly
with a small error (typically of size around 10-15 or 10-16). The analysis of the effect of round-
off error in calculations is a very interesting and technical subject, but it is not one of the things
that we shall do in this course
3. The quantities ex and log(x) (log to base e) are computed using the MATLAB commands exp
and log. Now compute e3 , log(e3) and eπ 163 .
4. Compute sin2 (π 6) + cos2 (π 6). (Note that typing sin^2(x) for sin2 x will produce an error!)
5. Use MATLAB as a calculator and try the following:
>>pi*pi – 10
>>sin(pi/4)
>> ans^2
>>xx=sin(pi/5 );
>>cos(pi/5 )
>>yy = sqrt(1 - xx*xx )
>>ans
6. Complex numbers are natural in MATLAB. Notice that the names of some basic operations
are unexpected, e.g., abs for magnitude. Try the following:
>>zz = 3 + 4i
>>conj(zz)

12 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Tutorial 2: Introduction Matrices and Graphics on MATLAB


Objective & Aim:
 To familiarize and practicing of the Matrices, Graphics on MATLAB software.

T2.1 Tutorial 2: Matrices, Graphics and File Creations


In this tutorial, you will learn
 how to create row and column vectors (one-dimensional arrays) and matrices (two-
dimensional arrays),
 how to do simple arithmetic operations (addition, subtraction, and multiplication by a
scalar) on arrays,
 how to do array operations, i.e. term-by-term multiplication (.*), division (./), and
exponentiation e.g. by a scalar (.^),
 how to use built-in functions with array arguments,
 how to set up and plot two-dimensional line plots,
 how to create m-files and execute it, and
 how to produce animated plots.

T2.1.1 Arrays
The basic data type in MATLAB is the rectangular array, which contains numbers (or even
strings of symbols). MATLAB allows the user to create quite advanced data structures, but for
many scientific applications it is sufficient to work with matrices.
These are two-dimensional arrays of data in the form
𝑎1,1 𝑎1,2 ⋯ 𝑎1,𝑚
A= 𝑎2,1 𝑎2,2 ⋯ 𝑎2,𝑚
𝑎𝑛,1 𝑎𝑛,2 ⋯ 𝑎𝑛,𝑚
Note that the notation implies that ai,j is the element of A which is stored in the ith row and jth
column. The number of rows n need not be the same as the number of columns m. We say that
such a matrix A is of dimension nxm. When n = 1 we have a column vector and when m = 1 we
have a row vector. When n = m = 1, A has a single entry, and is commonly known as a scalar.
You do not need any prior knowledge of matrices to follow the following part of notes. A
matrix can be created in MATLAB using any name starting with a letter and containing only
letters, numbers and underscores (_). A major distinction between MATLAB and many other
programming languages is that in MATLAB one does not have to declare in advance the size of
an array. Storage of arrays is dynamic in the sense that an array with a given name can change
size during a MATLAB session.
The MATLAB command

13 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

>>x = exp(1); creates a single number whose value is e and stores it in the variable whose name
is x. This is in fact a 1x1 array.
The command
>>y = [0,1,2,3,4,5];
creates a 1x6 array whose entries are equally spaced numbers between 0 and 5 and stores it in the
variable y. A quicker way of typing this is y = [0:5]. More generally, the command [a:b:c] gives
a row vector of numbers from a to c with spacing b. If the spacing b is required to be 1, then it
can be left out of the specification of the vector. Matrices can change dimension, for example if x
and y are created as above, then the command
>>y = [y,x]
Produces the response
y=
Columns 1 through 7
0 1.0000 2.0000 3.0000 4.0000 5.0000 2.7183
In the example above, we created a row vector by separating the numbers with commas. To
get a column vector of numbers you separate them by semi-colons, e.g. the command
>> z = [1;2;3;4;5]
z=
1
2
3
4
5
creates a 5x1 array of equally spaced numbers between 1 and 5.
>> size(z)
ans =
5 1

Matrices of higher dimension are created by analogous arguments. For example, the command

14 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

>> A = [[1,2,3];[4,5,6];[7,8,10]]
1 2 3
puts the matrix into the variable A. Similarly, [4 5 6]
7 8 9
>>C = [16 3 2 13; 5 10 11 8; 9 6 7 12; 4 15 14 1]
C=
16 3 2 13
5 10 11 8
9 6 7 12
4 15 14 1
The transpose of C can be obtained as
>> C'
ans =
16 5 9 4
3 10 6 15
2 11 7 14
13 8 12 1
The sum of the elements on the main diagonal is obtained with the sum and the diag functions.
sum(diag(C))
ans =
34
T2.1.2 Simple Operations
If A and B are two matrices (which must be of the same size)
>> A = [[1,2,3];[4,5,6];[7,8,10]];
>> B = [[11,12,13];[14,15,16];[17,18,19]];
then
>> A+B
ans =

15 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

12 14 16
18 20 22
24 26 29
A.*B produces a new matrix, the entries of which are the products of the corresponding entries
of A and B. That is,
A.*B
ans =
11 24 39
56 75 96
119 144 190
Note that the . before the * indicates that the multiplication is done pointwise (or elementwise).
Similarly the pointwise quotient is A./B. Hence it is called scalar operation.
Some of you may also know about the matrix product AB which is defined even when A and B
are not of the same size, as long as the number of columns of A is the same as the number of
rows of B. This (quite different) product is obtained in MATLAB with the command A*B. That
is,
>> A*B
ans =
90 96 102
216 231 246
359 384 409
Every function can be applied to a whole matrix (or more generally array) of numbers. The result
will be an array of the same size, each element of which is the result of the function applied to
the corresponding element of the array or matrix. For example, if
>>y=[0:5];
Then the command
>>z = log(y)
z=
Columns 1 through 6
-Inf 0 0.6931 1.0986 1.3863 1.6094

16 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Note that log is the built-in function which computes loge (or ln) and that the resulting vector y
above is the result of applying loge to each entry of the previous vector y. MATLAB handles an
infinite answer by recording it as Inf and producing a warning. Note that these calculations verify
that loge(1) = 0 and loge(e) = 1.
T2.1.3 MATLAB Build-in Functions with Array
Vectors and matrices can be created in your MATLAB workspace by various methods. For
example there are many built-in commands: zeros(n,m) and ones(n,m) produces nxm matrices of
zeros and ones.
For example:-
>> zeros(3,4)
ans =
0 0 0 0
0 0 0 0
0 0 0 0
The command rand(n,m) produces a random nxm matrix. The command eye(n,n) produces the
nxn identity matrix with ones on its diagonal (from top left to bottom right) and zeros elsewhere.
Try these out for some values of n and m.
For example;-
>> rand(3,2)
ans =
0.8147 0.9134
0.9058 0.6324
0.1270 0.0975
>> rand(3,2)
ans =
0.2785 0.9649
0.5469 0.1576
0.9575 0.9706
>> eye(4,4)
ans =

17 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
N.B:-Type help command for each of these commands and read the MATLAB manual page for
each of them.
If the entries of a vector or matrix are given by a formula then there is a simple way to create
them, called for loop. For example, if we wish to create a vector x whose entries are the squares
of the natural numbers from 1 to 10 we can use the command
for :-Repeat statements a specific number of times.
The general form of a for statement is:
for variable = expr, statement, ..., statement END
>> for i=1:10; x(i) = i^2; end;x
x=
Columns 1 through 10
1 4 9 16 25 36 49 64 81 100

T2.1.4 Graphics
The simplest plotting command is plot. If x and y are any two vectors of the same size, then the
command
>>x=[0:10];
>>y=[0:0.1:1];
>>plot(x,y)
>>plot(x,y,'r*')
>>stem(x,y)
The resultant plots of all three are shown in Figure T2.1, T2. 2, and T2. 3 which are the
components of y in turn against the components of x, using different mediums.

18 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Figure T2. 1 Graphical presentation using plot(x,y) function of the MATLAB

Figure T2. 2 Graphical presentation using plot(x,y,’r*’) function of the MATLAB

Figure T2. 3 Graphical presentation using stem(x,y) function of the MATLAB.

19 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

There are several variations of this command. Type help plot or doc plot --> Graphics -->
Line plotting to find out more about this command. Other simple plotting commands, for which
you should consult the help pages --> Graphics --> 2-D Plots, include fplot, bar, hist, stem, stairs
and pie.
The example below produces three simple plots. The first is a static plot comparing the
function sin(x) with its approximation x-x3/3. The second is a static plot producing a helix
(spiral) and the third plot is dynamically using the comet command. The plots come up in the
Figure window. They can be printed or copied and pasted into another document if required.
>> % program to explore the use of some simple plotting commands
>> % First example discusses a cubic approximation to sin(x)
>> x = pi*([0:100]/100); % row vector of 101 equally spaced points between 0 and pi
>> hold off % means that the next plot will be drawn on a clean set of
axes
>> plot(x,sin(x),'r*') % plots sin(x) against x using red *'s

Figure T2. 4 Graphical presentation by plot function of the MATLAB.


>> hold on % next plot will be drawn on top of the present one
>> plot(x,(x-x.^3/factorial(3)),'bo') % plots values of x - (x^3)/3! against x using blue o's

20 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Figure T2. 5 Graphical presentation by plot function of the MATLAB.


>> pause % causes execution to stop until you next hit
% any key (useful for admiring intermediate
plots)
% second example draws a 3D helix composed of the points
% (sin(t),cos(t),t) for t ranging over
% 301 equally spaced points between 0 and 100
>> t=[0:300]/3;
>> hold off
>> plot3(sin(t),cos(t),t,'r*',sin(t),cos(t),t,'g')
The result is shown in Fig. T2.6.

21 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Figure T2. 6 Graphical presentation by plot function of the MATLAB.

T2.1.5To Create m-files and Execute it


MATLAB is an interpretive language, i.e., commands typed at the MATLAB prompt are
interpreted within the scope of the current MATLAB session. However, it is tedious to type in
long sequences of commands each time MATLAB is used to perform a task.
There are two means of extending MATLAB’s power: scripts and functions. Both make use of
m-files (named because they have a .m extension and therefore they are also called dot-m files)
created with a text editor. The advantage of m-files is that commands are saved and can be easily
modified without retyping the entire list of commands.
Here we make use of only the MATLAB script files, which are sequences of commands typed
with an editor and saved in an m-file. To create an m-file, open the Edit Window and type the
program. One example of this is shown in Figure T2. 1. Top of the Edit Window shows the file
name with its path. By default, its path is set as D:\dsp lab\mat lab exp\tutoria1212.m, where
tutoria1212.mis the name of m-file.

22 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Figure T2. 7 Edit Window with a program to demonstrate the animated plot
The Edit Window is essentially a programming text editor, with the features of the languages
of MATLAB being highlighted in different colors. Comments in m-files appear in green,
variables and numbers appear in black, character strings appear in red, and language keywords
appear in blue. Do not create a m-file with the same name as a MATLAB function or command.
To execute the m-file, either type the file name in the command prompt (>>) or press F5. The m-
file tutoria1212.m gives the figure which is shown in Figure T2. 2.
A second type of m-file is a function file which is generated with an editor exactly as the script
file. Students are advised to learn this too.

Figure T2. 8 Animated plot of m-file tutoria1212.m

T2.1.6 Animated Graphs


We have seen an animated plot in Figure T2.9. It plots the waves given by the functions
sin(x − t), cos(x − t) for values of x between 0 and 8π and for time t = 0; 1; 2; ….; 100. Here afor

23 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

loop is used to advance time in an obvious way. Note the very important draw now command,
which causes the plot to be put on the screen at each time step, thus allowing the animation to be
seen.
Let us see some more animated plot of set of values using the command ``comet''
>> t = [0:300]/3;
>> comet3(sin(t),cos(t),t)
The resultant plot is shown in Figure T2. 10

Figure T2. 9 Animated plot using command ‘comet3’

Exercises for Tutorial 2 -


1. Create 3x3 matrix whose entries in the first row and first column are all equal to 1 and all the
other entries are equal to e.
2. Plot the function exp(-x) sin(8x) for 0 ≤ x ≤ 2π .
3. Generate a row vector x of even integers, each component x(i) of which lies in the range 0 ≤
x(i) ≤10 . Compare the effects of the operations
(i) sin(x), (ii) x.^2, (iii) x.^(-1).
4. Generate the 1x10 row vector v whose ith component is cos(iπ/4).
5. Modify the script file ani.m in the editor/debugger window to do an animated plot of the
functions sin(4(x-t)) and 4 cos(x-t) for x between 0 and 8π and t from 0 to 10. Experiment with
different line colours and types.

24 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

6. Plotting is easy in MATLAB, for both real and complex numbers. The basic plot command
will plot a vector yy versus a vector xx. Try the following:
>>xx = [-3 -1 0 1 3 ];
>>yy = xx.*xx - 3*xx;
>>plot( xx, yy )
>>zz =xx + yy*sqrt(-1);
>>plot( zz )
Drop the semicolons, if you want to display the values in the xx, yy, and zz vectors. Use help to
learn how the operation xx.*xx works; compare to matrix multiply.
7. Be sure to understand the colon notation. In particular, observe what the following MATLAB
code will produce:
>>jk1= 2 : 4 : 17
>>jk1= 99 : -1 : 88
>>ttt= 2 : (1/9) : 4
>>tpi=pi * [2 : (-1/9) : 0 ]
8. Extracting or inserting numbers in a vector is very easy to do. Consider the following
definition:
>>xx = [ ones(1,4), [2:2:11], zeros(1,3) ]
>>xx(3:7)
>>length(xx)
>>xx(2:2:length(xx))
Analyse the result echoed from the last three lines of this code.

25 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Experiment Part one: Generation of Discrete-Time Signals in Time Domain


Lab #1: Generation of Discrete-Time Signals in Time Domain
1. Objectives:
 To generate unit step, impulse, exponential and sinusoidal discrete-time signals.
 To study the characteristics of these signals.

2. Materials Required: Matlab Software and Computers.


3. Background Theory
In this part of the experiment, some discrete-time signals will be generated and studied. These
signals may be applied in real-time digital signal processing.
(a) The impulse discrete signal
The discrete-time impulse signal can be represented as
𝛿 (𝑛) = 1, 𝑛 = 0
0, 𝑛 ≠ 0

b) The unit step discrete signal

26 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Mathematically, a unit step discrete signal is written as:


𝑢(𝑛) = 1, 𝑛 ≥ 0
0, 𝑛 < 0
(c) The ramp discrete signal
Mathematically, the ramp discrete signal is written as
𝑟(𝑛) = 𝑛, 𝑛 ≥ 0
0, 𝑛 < 0

(d) The exponential discrete signal


Mathematically, the exponential discrete signal is written as

𝑥(𝑛) = 𝐴𝛼 𝑛
Α is real value. 𝛼 Can be real or complex.
If 0≤ 𝛼 ≤ 1, then the signal x (n) will decay exponentially. If 𝛼 ≥ 1, then the signal x (n) will
grow without bound.
(e) The sinusoidal discrete signal
Mathematically, the sinusoidal discrete signal is written as
𝑥(𝑛) = 𝐴𝑐𝑜𝑠(𝑤𝑜 𝑛 + 𝜑)
Where A is the amplitude, 𝑤𝑜 is the angular frequency and 𝜑 is the phase.
(f) RANDOM SIGNALS
A random signal of length N with samples uniformly distributed in the interval [0,1] can be
generated using the following MATLAB command
𝑥 = 𝑟𝑎𝑛𝑑(1, 𝑁)
Similarly a random signal of length N with samples normally distributed with zero mean and
variance one can be generated using the MATLAB command.
𝑥 = 𝑟𝑎𝑛𝑑𝑛(1, 𝑁)

4. Procedures:

27 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

The procedures to develop a Matlab code to generate different discrete-time signals are stated
below for each signal type separately.

(a) Impulse response


 Define the starting and ending points on the graph to plot impulse response.
 Enter the starting and ending points.
 Enter the location of the impulse response.
 Plot the impulse response.
The following program 1.1 is written to simulate the unit sample response function. Write the
code on Matlab M-file and run it for different input values

Code 1:-
% program 1.1 to simulate the
impulse/unit sample response function
n = -10:20;
n0=0;
% Generate the unit sample sequence
u =n-n0==0;
% Plot the unit sample sequence
stem(n,u);
xlabel('Time index n');
ylabel('Amplitude');
title('Unit Sample Sequence');
axis([-10 20 0 2]);

Figure P1L1. 1 Impulse response-sample output of program1.1

28 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Code#2:-
% program 1.2 to simulate the impulse/unit sample
response function by using input keyword
clc;
clf;
N=input('unit impulse sequence of length N=');
M=input('unit impulse sequence delayed by M=');
n=-N:M;
% Generate the unit sample
% sequence
u=[zeros(1,N) 1 zeros(1,M)];
% plot of unit sample sequence
stem(n,u)
xlabel('Time index n');
ylabel('Amplitude');
title('unit Sample Sequence');
axis([-10 20 0 1.2]);
grid on;

N.B:- input prompt for user input.


RESULT = input(PROMPT) displays the PROMPT string on the screen, waits for input
from the keyboard, evaluates any expressions in the input, and returns the value in RESULT.
Take the following sample inputs to get the display of figure 1.2:
Unit impulse sequence of length N=10 Unit impulse sequence delayed by M=5

Figure P1L1. 2 Impulse response-sample output of program1.2

(b) Step response

29 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

 Define the interval of the unit step function, that is, enter the starting and ending points.
 Enter the shift point no.
 Plot the sequence.
Write the following code on Matlab M-file and run it.
Code#2:-
% Program 1.3 to generate Unit Step Signal
clc; clear all; close all hidden;
N =input('Enter the length of the sequence (even number) N
: ');
m =input('Enter the location of impulse sample
m : ');
n=[(-N/2):(N/2)];
unit_step=[n>=m];
stem(n,unit_step);
xlabel('Time index n');
ylabel('Amplitude');
title('Unit Step Response')

Enter the length of the sequence (even number) N : 10


Enter the location of impulse sample m: 2

Figure P1L1. 3 The unit step response: sample output of program 1.3

C) Unit Ramp Signal


Step 1: Start the program.
Step 2: Enter the length of the ramp sequence.
Step 3: Plot the sequence in discrete form.

30 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Step 4: Stop the program.


Program for Unit Ramp Signal
%Program P1.4
%can be used to generate and plot a unit step sequence
clc;
clf;
N=input('unit impulse sequence of length N=');
M=input('unit impulse sequence delayed by M=');
n=-N:M;
% Generate the unit sample
% sequence
u=[zeros(1,N) 1 ones(1,M)];
% plot of unit sample sequence
stem(n,u)
xlabel('Time index n');
ylabel('Amplitude');
title('unit Sample Sequence');
axis([-10 20 0 1.2]);
grid on;

unit ramp sequence of length N=10


unit ramp sequence delayed by M=20

Figure P1L1. 4 Discrete-time representation of the unit ramp signal.

(d) Exponential Sequences


 Define the exponential function (𝑥(𝑛) = 𝐴𝛼 𝑛 )
 Fix the values of its amplitude A, base , and time index n as shown in program 1.4.
 Write a matlab program 1.4 and run it.
 Observe the plots.

31 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

% Program 1.5
% Generation of a complex exponential sequence
clf;
c = -(1/12)+(pi/6)*i; % define the base which is complex
A = 2; % amplitude of the function
n = 0:40; % time index range
x = A*exp(c*n); % the function defined
subplot(2,1,1);
stem(n,real(x)); % plot real part of x
xlabel('Time index n');ylabel('Amplitude');
title('Real part');
subplot(2,1,2);
stem(n,imag(x)); % plot imaginary part
xlabel('Time index n');ylabel('Amplitude');
title('Imaginary part');

N.B:- subplot Create axes in tiled positions.


H = subplot(m,n,p), or subplot(mnp), breaks the Figure window into an m-by-n matrix of small
axes, selects the p-th axes for the current plot, and returns the axes handle.

Figure P1L1. 5 Complex exponential sequence: sample plot of program 1.5.


(d) Sinusoidal sequences
 Define the sinusoidal sequence function (𝑥(𝑛) = 𝐴𝑐𝑜𝑠(𝑤𝑜 𝑛 + 𝜑)) with its
amplitude

32 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

 Angular frequency 𝑤𝑜 , phase constant 𝜑 .


 Write the program 1.6 in matlab.
 Observe the outputs/plots of the program.
% Program 1.5
% Generation of a sinusoidal sequence
n = 0:40;
f = 0.1;
phase = 0;
A = 1.5;
arg = 2*pi*f*n - phase;
x = A*cos(arg);
clf; % Clear old graph
stem(n,x); % Plot the generated sequence
axis([0 40 -2 2]);
grid;
title('Sinusoidal Sequence');
xlabel('Time index n');
ylabel('Amplitude');
axis;

Figure P1L1. 6 Sinusoidal sequence: sample plot of program 1.6.

e) Random Signals
One of the applications of digital signal processing is to remove random noise from a
signal which has been corrupted by noise. Let s [n] be the signal corrupted by a random noise d
[n] resulting in the noisy signal x[n] = s[n] + d[n]. The objective is to operate on x [n] to
generate a signal y [n] which is a reasonable approximation to s [n]. This can be done by
averaging the samples of x[n]. This is referred to as moving average filter
𝑦[𝑛] = 1/3(𝑥(𝑛 − 1) + 𝑥[𝑛] + 𝑥(𝑛+))

33 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

%Program P1.7
%signal Smoothing by Averaging
clc;
clf;
R=51;
d=0.8*(rand(R,1)-0.5);
m=0:R-1;
s=2*m.*(0.9.^m);
x=s+d';
subplot(2,1,1);
plot(m,d,'r',m,s,'g--',m,x,'b-.'); grid on;
xlabel('Time index n');ylabel('Amplitude');
legend('d[n] ','s[n] ','x[n] ');
x1=[0 0 x];x2=[0 x 0];x3=[x 0 0];
y=(x1+x2+x3)/3;
subplot(2,1,2);
plot(m,y(2:R+1),'r-',m,s,'g--');grid on
legend('y[n] ','s[n] ');
xlabel('Time index n');ylabel('Amplitude');

Figure P1L1. 7 Sinusoidal sequence: sample plot of program 1.7

5. Lab Task
(i) Impulse and unit-step responses
1. Run Programs 1.1 and 1.3 to generate the unit sample and unit step sequences for
different input parameters.
2. What are the purposes of the commands input, stem, and title?
3. Write a matlab code that displays graphically the following discrete-time signals:
(a) x (n)=u(n+4)-u(n-4)

34 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

(b) x(n)=δ(n)+δ(n+2)
(ii) Exponential sequences
1. Run Program 1.4 and generate the complex-valued exponential sequence.
2. Which parameter controls the rate of growth or decay of this sequence? Which
parameter controls the amplitude of this sequence?
3. What are the purposes of the operators’ real and imag?
4. What is the purpose of the command subplot?

(iii). Write MATLAB programs to generate the square, triangular and sawtooth wave sequences
of the types shown in Figures below.
1. square wave

Figure P1L1. 8 square wave


2. Sawtooth Wave Sequence

35 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Figure P1L1. 9 Sawtooth Wave Sequence

Experiment Part Two: Simulation of Discrete-Time Systems


Lab#1The Moving Average System
1. Objectives:
 To study the noise filtering characteristics of the moving average system.

2. Background Theory
A causal version of the three-point smoothing filter is given by
𝑦[𝑛] = 1/3(𝑥[𝑛] + 𝑥 [𝑛 − 1] + 𝑥[𝑛 − 2])
Generalizing the above equation we obtain
𝑀−1
1
𝑦[𝑛] = ∑ 𝑥[𝑛 − 𝑘]
𝑀
𝑘=0

36 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

This defines a causal M-point smoothing FIR filter. The above equation is also known as a
moving average filter. Program 2.1 illustrates its use in filtering high-frequency components
from a signal composed of a sum of several sinusoidal signals.

3. Materials Required: Matlab Software and Computers.


4. Procedure:
 Define the time index n (example 0 to 100).
 Define a low frequency sinusoidal signal s1 and high frequency component signal s2.
𝑠1 = cos(𝑤1 𝑡) and 𝑠2 = cos(𝑤2 𝑡); 𝑤1 𝑎𝑛𝑑 𝑤2 are angular frequencies ; let
𝑤1 = 0.05, 𝑤2 = 0.47
 Add the two signal frequencies: x=s1+s2.
 Input the desired length of the filter.
 Plot each signals s1, s2 and added signal x.
 Filter out the noisy signal x and implement this action using the filter command as shown
in program 2.1.
 Observe the outputs from the plots.

% Program 2.1
% Simulation of an M-point Moving Average Filter
% Generate the input signal
n = 0:100;
s1 = cos(2*pi*0.05*n); % A low frequency sinusoid
s2 = cos(2*pi*0.47*n); % A high frequency sinusoid
x = s1+s2;
% Implementation of the moving average filter
M = input('Desired length of the filter = ');
num = ones(1,M);
y = filter(num,1,x)/M;
% Display the input and output signals
clf;
subplot(2,2,1);
plot(n,s1);
axis([0, 100, -2, 2]);
xlabel('Time index n'); ylabel('Amplitude'); title('Signal # 1');
subplot(2,2,2);
plot(n,s2);
axis([0, 100, -2, 2]);
xlabel('Time index n'); ylabel('Amplitude'); title('Signal # 2');
subplot(2,2,3);
plot(n,x);
axis([0, 100, -2, 2]);
xlabel('Time index n'); ylabel('Amplitude'); title('Input Signal');
subplot(2,2,4);
plot(n,y);
axis([0, 100, -2, 2]);
xlabel('Time index n'); ylabel('Amplitude'); title('Output Signal');
axis;

37 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Figure P2L1. 1 The illustration of the moving average filter Program 2.1 .

5. Lab Task:-
1. Run the above program for M = 2 to generate the output signal with x[n] = s1[n]+ s2[n] as the
input. Which component of the input x[n] is suppressed by the discrete time system simulated by
this program?
2. If the system is changed from y[n] =0.5(x[n] + x[n - 1]) to y[n] =0.5(x[n] - x[n - 1]), what
would be its effect on the input x[n] = s1[n] + s2[n]?
3. Run Program 2 .1 for other values of filter length M, and various values of the frequencies of
the sinusoidal signals s1 [n] and s2 [n]. Comment on your results.

38 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Lab #2 Discrete-Time Systems in the Time Domain for LTI


1. Objective:
 To simulate and study the characteristics of the linearity property of discrete-time
systems.
 To investigate the time-invariance property of a discrete-time system
2. Background theory:
Mathematically, a discrete-time system is described as an operator T[.] that takes a sequence x(n)
called excitation and transforms it into another sequence y(n) (called response). Discrete time
systems can be classified into two categories i) LTI systems ii) NON-LTI systems. A discrete
system T[.] is a linear operator T[.] if and only if T[.] satisfies the principle of superposition,
namely

𝑇[𝑎1 𝑥1 (𝑛) + 𝑎2 𝑥2 (𝑛)] = 𝑎1 𝑇[𝑥1 (𝑛)] + 𝑎1 𝑇[𝑥2 (𝑛)]

A discrete system is time-invariant if shifting the input only causes the same shift in the output.

A system is said to be bounded-input bounded-output (BIBO) stable if every bounded input


produces a bounded output.
|𝑥(𝑛)| < ∞ → |𝑦(𝑛)| < ∞, ∀𝑥, 𝑦

An LTI system is BIBO stable if and only if its impulse response is absolutely sum able.

39 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

𝐵𝐼𝐵𝑂 𝑆𝑡𝑎𝑏𝑖𝑙𝑖𝑡𝑦 ↔ ∑|ℎ(𝑛)| < ∞


A system is said to be causal if the output at any instant depends only on the present & past
values only.

An LTI system is causal if and only if the impulse response is

ℎ(𝑛) = 0, 𝑛 < 0
3. Materials Required: Matlab Software and Computers.

4. Procedure:
a. Linearity and Non-Linearity:
We now investigate the linearity property of a causal system of described by the following
equation.
Ex:-
𝑦[𝑛] − 0.4𝑦[𝑛 − 1] + 0.75𝑦[𝑛 − 2] = 2.2𝑥[𝑛] + 2.3𝑥[𝑛 − 1] + 2.4𝑥[𝑛 − 2]

%Program P2.2
clc
n = 0:40;
a = 2; b = -3;
x1 = cos(2*pi*0.1*n);
x2 = cos(2*pi*0.4*n);
x = a*x1 + b*x2;
num =[1 0.3 2.4] ;
den = [1 0 -1 0.75 3];
ic = [0 0 0 0]; % Set zero initial conditions
y1 = filter(num,den,x1,ic); % Compute the output y1[n]
y2 = filter(num,den,x2,ic); % Compute the output y2[n]
y = filter(num,den,x,ic); % Compute the output y[n]
yt = a*y1 + b*y2;
d = y - yt; % Compute the difference output d[n]
% Plot the outputs and the difference signal
subplot (4,1,1)
stem(n ,y);
ylabel('Amplitude');title('Output Due to Weighted Input');
subplot(4,1,2)
stem(n,yt);
ylabel('Amplitude');title('Weighted Output');
subplot(4,1,3)
stem(n,d);
xlabel('Time index n');ylabel('Amplitude');title('Difference Signal');
subplot(4,1,4)
plot(n,d);
xlabel('Time index n');ylabel('Amplitude');title('Difference Signal');

Following program simulates the above mentioned equation.

40 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Figure P2L1. 2 Demonstration of linearity property of the discrete-time system .

b. Time-Invariant and Time-Varying Systems:


We next investigate the time-invariance property. Following program simulates following
difference equation
𝑦[𝑛] − 0.4𝑦[𝑛 − 1] + 0.75𝑦[𝑛 − 2] = 2.2𝑥[𝑛] + 2.3𝑥[𝑛 − 1] + 2.4𝑥[𝑛 − 2]

%Program P2.3
clc
n = 0:40; D = 10;a = 3.0;b = -2;
x = a*cos(2*pi*0.1*n) + b*cos(2*pi*0.4*n);
xd = [zeros(1,D) x];
num = [2.2 2.3 2.4];
den = [1 -0.4 0.75];
ic = [0 0];% Set initial conditions
% Compute the output y[n]
y = filter(num,den,x,ic);
% Compute the output yd[n]
yd = filter(num,den,xd,ic);
% Compute the difference output d[n]
d = y - yd(1+D:41+D);
% Plot the outputs
subplot(3,1,1)
stem(n,y);
ylabel('mplitude');title('Output y[n]');grid;
subplot(3,1,2)
stem(n,yd(1:41));
ylabel('Amplitude');title(['Output Due to Delayed Input x[n ,
num2str(D),]']);grid;
subplot(3,1,3)
stem(n,d);
xlabel('Time index n'); ylabel('Amplitude');title('Difference
41Signal');grid;
| P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN
Digital Signal Processing Lab Manual

Figure P2L1. 3 Demonstration of Time-Invariant and Time-Varying Systems of the discrete-time


system

C.Stability of LTI Systems


As indicated, an LTI discrete-time system is BIBO stable if its impulse response is absolutely
sum able. It therefore follows that a necessary condition for an IIR LTI system to be stable is
that its impulse response decays to zero as the sample index gets larger. Program P4 8 is a
simple MATLAB program used to compute the sum of the absolute values of the impulse
response samples of a causal IIR LTI system. It computes N samples of the impulse response
sequence, evaluates
𝑘

𝑠(𝑘) = ∑ |ℎ[𝑛]|
𝑛=0
for increasing values of K, and checks the value of |h[K]| at each iteration step. If the value of
|h[K]| is smaller than 10-6, then it is assumed that the sum S(K) has converged and is very
close to S(∞).

42 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

%program p2.3
clc;clf;
num=[1 -0.8]; den=[1 1.5 0.9];
N=200;
h=impz(num,den,N+1);
parsum=0;
for k=1:N+1
parsum=parsum+abs(h(k));
if abs(h(k) <10^(-6)),break, end
end
n=0:N;
stem(n,h);
xlabel('Time index n');ylabel('Amplitude');
disp('Value=');
disp(abs(h(k)));

Figure P2L1. 4 Demonstration Stability of LTI Systems of the discrete-time system.


5. Lab Task:-
1. Run Program 2.3 and compare the output sequences y[n] and yd[n - 10]. What is the
relation between these two sequences? Is this system time-invariant?
2. Modify given program to compute the output sequences y1[n], y2[n], and y[n] of the
above system. Compare y[n] with y[n]. Are these two sequences equal? Is this system
linear?
A.𝑦[𝑛] − 𝑦[𝑛 − 2] + 0.75𝑦[𝑛 − 4] = 2.2𝑥[𝑛] + 2.3𝑥[𝑛 − 3] + 2.4𝑥[𝑛 − 5]
B. 𝑦[𝑛] − 0𝑦[𝑛 − 2] + 1𝑦[𝑛 − 4] = 𝑥[𝑛] + 0𝑥[𝑛 − 3] + 𝑥[𝑛 − 4]
C.𝑦[𝑛] − 𝑦[𝑛 − 2] + 0.75𝑦[𝑛 − 4] + 3𝑦[𝑛 − 5] = 𝑥[𝑛] + 0.3𝑥[𝑛 − 1] + 2.4𝑥[𝑛 − 3]
3. Modify program to simulate the above system and determine whether this system is time-
invariant or not.
𝐴. 𝑦[𝑛] − 𝑦[𝑛 − 1] + 𝑦[𝑛 − 2] + 0.01𝑦[𝑛 − 4] = 𝑥[𝑛] + 2.3𝑥 [𝑛 − 1] + 2.4𝑥[𝑛 − 5]

43 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

𝐵. 𝑦[𝑛] − 𝑦[𝑛 − 2] = 𝑥[𝑛] + 0.02𝑥[𝑛 − 2] + 𝑥 [𝑛 − 4]


𝐶. 𝑦[𝑛] = 𝑥[𝑛] + 2.3𝑥[𝑛 − 3] + 2.4𝑥[𝑛 − 5]
4. Modify program to simulate the above system and determine whether this system is
Stability of LTI Systems or not.
𝐴. 𝑦[𝑛] = 𝑥(𝑛 − 1) + 1/2(𝑥[𝑛 − 2]) + 𝑥[𝑛]

Lab #3 Linear Convolution


1. Objective:
To compute the linear convolution of a linear time invariant (LTI) system from it’s a finite length
Impulse response h (n) and finite length sequence x (n).
2. Background Theory:
The output y (n) to the input x(n) of the LTI system with impulse response h(n) is given by

𝑦(𝑛) = ℎ(𝑛) ∗ 𝑥(𝑛) = ∑ 𝑥(𝑘)ℎ(𝑛 − 𝑘)


𝑛=−∞

This equation is termed as the convolution sum.


To find the output of a discrete system y(n) to an input x(n), we need the impulse response,
h(n),
For the system. h(n) is the output of the system when the input is δ(n), where δ(n) is the impulse
Signal.
The convolution operation is implemented in MATLAB by the command conv, provided the
two
Sequences to be convolved are of finite length. For example, the output sequence of a finite
impulse response system can be computed by convolving its impulse response with a given
finite-length input sequence. The following MATLAB program illustrates this approach.

3. Materials Required: Matlab Software and Computers.


4. Procedure:
 Take the input x[n] and the impulse response h[n] of the linear time-invariant system.
 Use the Matlab function conv to produce the convolved output of the system.
 Use the filter function to produce the output.

44 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

 See the difference between the one obtained by filter function and the conv function.

% Program P1.1
clf;
h = [3 2 1 -2 1 0 -4 0 3]; % impulse response
x = [1 -2 3 -4 3 2 1]; % input sequence
y = conv(h,x);
n = 0:14;
subplot (2,1,1);
stem(n,y);
xlabel('Time index n'); ylabel('Amplitude');
title('Output Obtained by Convolution');grid;
x1 = [x zeros(1,8)];
y1 = filter(h,1,x1);
subplot(2,1,2);
stem(n,y1);
xlabel('Time index n'); ylabel('Amplitude');
title('Output Generated by Filtering');grid;

Figure P2L1. 5 linearly convolved output of LTI system to a finite length input.

5. Lab Task
1. Run Program 1.1 to generate y[n] obtained by the convolution of the sequences h[n] and x[n],
and to generate y1[n] obtained by filtering the input x[n] by the FIR filter h[n]. Is there any
difference between y[n] and y1[n]? What is the reason for using x1[n] obtained by zero-padding
x[n] as the input for generating y1[n]?
2. Modify Program 1.1 to develop the convolution of a length-15 sequence h[n] with a length-10
sequence x[n], and repeat question 1. Use your own sample values for h[n] and x[n].

45 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Lab #3:- Circular Convolution.


1. Objective
 To write a program in MATLAB to perform the circular convolution for testing the
periodic systems.

2. Background Theory
The linear convolution possesses difficulty because of the finite word length of the real-time
processors. This difficulty has been solved by the introduction of the circular convolution. The
circular convolution is otherwise called the periodic convolution and it gets repeated at every
sample period, N. Hence, by definition the periodic convolution is written as

𝑦(𝑛) = 𝑥(𝑛) ∗ ℎ(𝑛) = ∑ 𝑥(𝑙 )ℎ(𝑛 − 𝑙)𝑁


𝑙=−∞
The suffix N indicates that it is periodic with N number of samples. The length of the circular
convoluted sequence will be the same as that of the length of the given sequences, whichever has
the higher length. In other words, the length of the convoluted sum will be either M or N
whichever is higher, where M be the length of the input sequence and N be the length of the
impulse response.

3. Materials Required: Matlab Software and Computers.


4. Procedure:
Step 1: Start the program.
Step 2: Enter the input sequence.
Step 3: Enter the impulse response.
Step 3: Get the convoluted output and plot the sequence in discrete form.
Step 4: Stop the program.

46 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Program for Circular Convolution


Input Condition:
Enter the input sequence x(n): [1 2 3 5 6]
Enter the impulse response h(n): [1 2 5 8]

% Program to perform Circular Convolution %


clc; close all hidden; clear all;
x=input('Enter the input seqence x(n): ');
h=input('Enter the impulse response h(n): ');
N1=length(x);
N2=length(h);
N=max(N1,N2);
N3=N1-N2;
% loop for getting equal length sequence
if(N3==0)
x=[x,zeros(1,N3)];
h=[h,zeros(1,N3)];
end
if(N1>N2)
h=[h,zeros(1,N3)];
end
if(N1<N2)
x=[x,zeros(1,-N3)];
end
% computation of circular convolved sequence
for n=1:N,
y(n)=0;
for i=1:N,
j=n-i+1;
if(j==0)
j=N+j;
end
if(j<0)
j=N+j;
end
y(n)=y(n)+x(i).*h(j);
end
end
display('Cicular Convoluted sum: ');
y
p=0:(N-1);
stem(p,y), grid;

47 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Figure P2L1. 6 Convoluted sum by using the circular convolution.

5. Lab Task:-
1. Perform the circular convolution of the following. Consider the first element of each sequence
at the index value zero.

A, x(n)=[2,1,2,1] and h(n)=[1,2,3,4]


B, x(n) =cos(2*pi*13n) where n=1:5 and h(n)= [1,2,3,4]
2. Modify Program 1.1 to develop the convolution of a length-15 sequence h[n] with a length-10
sequence x[n], and repeat question 1. Use your own sample values for h[n] and x[n].

48 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Experiment Part Three: The Frequency Domain Analysis


In this experiment part, the frequency domain analysis of signals will be carried out by using
different analysis methods such as discrete time Z- transform Fourier transform, discrete Fourier
transform, and fast Fourier transform.

Lab #1: Discrete-Time Fourier Transform (DTFT) Computation


1. Objective
- To transform time domain signals to frequency domain by DTFT.

2. Back ground Theory


The discrete-time Fourier transform (DTFT) X(ejω) of a sequence x[n] is defined by

𝑗𝜔𝑡 )
𝑥(𝑒 = ∑ 𝑥[𝑛]𝑒 −𝑗𝜔𝑛
𝑛=−∞

Some class of DTFTs that can be expressed in rational form of equation xxx can be evaluated
using Matlab command freqz.
∑𝑀
𝑘=0 𝑏1 𝑒
−𝑗𝜔𝑡
𝑏0 + 𝑏1 𝑒 −𝑗𝜔𝑡 + ⋯ … . 𝑏𝑀 𝑒 −𝑗𝜔𝑡
𝑥(𝑒 𝑗𝜔𝑡 ) = =
∑𝑀
𝑘=0 𝑎𝑖 𝑒
−𝑗𝑘𝑡 𝑎0 + 𝑎1 𝑒 −𝑗𝜔𝑡 + ⋯ … . 𝑎𝑁 𝑒 −𝑗𝜔𝑡
Where k=0,1,2…..M, i=1,2,3…….N and ai and bk are constant coefficients.
In general X(ejω) is a complex function of the real variable ω and can be written as

𝑥(𝑒 𝑗𝜔𝑡 ) = 𝑥𝑟𝑒 (𝑒 𝑗𝜔 ) + 𝑗𝑥𝑖𝑚 (𝑒 𝑗𝜔 )


where Xre(ejω) and Xim(ejω) are, respectively, the real and imaginary parts of X(ejω), and are real
functions of ω. X(ejω) can alternately be expressed in the form

𝑥(𝑒 𝑗𝜔𝑡 ) = |𝑥 (𝑒 𝑗𝜔 )|𝑒 𝑗𝜃𝜔

where

𝜃(𝜔) = 𝑎𝑟𝑔{ 𝑥(𝑒 𝑗𝜔𝑡 )}

The quantity |X(ejω)| is called the magnitude function and the quantity θ(ω) is called the phase
function , with both functions again being real functions of ω. In many applications, the Fourier
transform is called the Fourier spectrum and, likewise, |X(e jω)| and θ(ω) are referred to as the
magnitude spectrum and phase spectrum, respectively.
The matlab command freqz gives the frequency response of a linear discrete-time system
specified by its transfer function.

49 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Input parameters:
 a vector b containing the coefficients of the numerator polynomial bl , for l =0, 1, . . .,M;
 a vector a containing the coefficients of the denominator polynomial , ai, for i =0, 1, . . .,
N;
 the number n of points around the upper half of the unit circle, or, alternatively, the
vector w of points in which the response is evaluated.
Output parameters:
 a vector containing the frequency response;
 a vector w containing the frequency points at which the frequency response is evaluated;
 With no arguments, freqz plots the magnitude and unwrapped phase of the transfer
function.

3. Materials Required: Matlab Software and Computers.


4. Procedure:
 Specify the frequency samples, vector w.
 Compute the DTFT of x(n) using equation (3.1a) and express in the form of equation
(3.1b).
 Obtain the numerator and denominator vectors: num and den.
 Use the matlab command freqz to return response with its real part, complex part and
phase response components.
 Plot the transformed signal using different matlab plotting commands.
 Use the developed program 3.1 and answer all the questions given next to the sample
outputs.

50 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

% Program 1.1
% Evaluation of the DTFT
clf;
% Compute the frequency samples of the DTFT
w = -4*pi:8*pi/511:4*pi;
num = [2 1];den = [1 -0.6];
H= freqz(num, den, w);
% Plot the DTFT
subplot(2,1,1)
plot(w/pi,real(H));grid
title('Real part of H(e^{j\omega})');xlabel('\omega /\pi');
ylabel('Amplitude');
subplot(2,1,2)
plot(w/pi,imag(H));grid
title('Imaginary part of H(e^{j\omega})');xlabel('\omega /\pi');
ylabel('Amplitude');
pause
subplot(2,1,1)
plot(w/pi,abs(H));grid
title('Magnitude Spectrum |H(e^{j\omega})|');xlabel('\omega /\pi');
ylabel('Amplitude');
subplot(2,1,2)
plot(w/pi,angle(H));grid
title('Phase Spectrum arg[H(e^{j\omega})]');xlabel('\omega /\pi');
ylabel('Phase, radians');

Figure P3: 1 Compute the frequency samples of the DTFT

51 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

5. Lab Task:-
1. What is the function of pause keyword?
2. Find the frequency response of the system with h (n) =n (0.1)n sin(n) for 0 ≤ ≤ 100. Plot
the magnitude and phase response of the system. Implement it in Matlab.
3. Find the Fourier transform of x(n) = δ(n) + (δ n – 1) + (n δ – 2) + (δ n – 3). Plot the
magnitude and the phase of X(jω).
4. Consider the difference equation 𝑦(𝑛) = 0.1𝑦(𝑛 − 1) + 0.2𝑦(𝑛 − 2) = 𝑥(𝑛) Find the
frequency response and plot its magnitude and phase.

52 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Lab #2:-The Discrete Fourier Transform (DFT) and Inverse DFT Computations
1. Objectives:
 To sample signals in frequency domain.
 To compute DFT of signals efficiently using fast Fourier Transforms (FFT) algorithms
embedded in Matlab.
 To compute circular convolution of two signals in frequency domain.
2. Background Theory
The N-point discrete Fourier transform (DFT) of a finite-length sequence x[n], defined for 0 ≤ n
≤ N − 1, is given by
𝑁−1

𝑋(𝐾 ) = ∑ 𝑥(𝑛)𝑊𝑁𝑘𝑛 , 𝑘 = 0,1,2 … . 𝑁 − 1


𝑛=0
Where 𝑊𝑁𝑘𝑛 −𝑗2𝜋/𝑁
=𝑒
The N-point DFT X[k] of a length-N sequence x[n], n = 0, 1. . . N−1, is simply the frequency
samples of its DTFT X(ejω) evaluated at N uniformly spaced frequency points, ω = ωk = 2πk/N,
k = 0, 1. . . N − 1, that is,
𝑋[𝐾 ] = 𝑋(𝑒 𝑗𝜔 )|𝑤=2𝜋𝑘/𝑁 k=0,1…….N-1
The N-point circular convolution of two length-N sequences g[n] and h[n], 0 ≤n ≤ N − 1, is
defined by
𝑀−1

𝑦(𝑛) = ∑ 𝑔[𝑚]ℎ[(𝑛 − 𝑚)]𝑁


𝑚=0
Where h [(n)] N = n modulo N.
Two important concepts used in the application of the DFT are the circular-shift of a sequence
and the circular convolution of two sequences of the same length.

The Fast Fourier Transform (FFT):


The Fast Fourier Transform, or the FFT, as it is popularly termed, is probably the single most
famous computer program in the field of electrical engineering. It is, essentially, a much faster
computation method of the DFT.
(a) Program to compute N-point DFT of N-point data sequence using direct method

3. Materials Required: Matlab Software and Computers.


4. Procedure:
 Specify the n-point data sequence x(n).
 Compute recursively the output X(k) as shown in program p1.1.
 Plot the magnitude of the computed sequence and observe its results.

53 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

a) perform Discrete Fourier Transform


% program 2.2
% Program to perform Discrete Fourier Transform %
clc; clear all; close all hidden;
x=input('The given sequence is x(n): ');
N=length(x);
for k=1:N
X(k)=0;
for n=1:N
X(k)=X(k)+x(n).*exp(-j.*2.*pi.*(n-1).*(k-1)./N);
end
X(k)
end
display('The DFT of the given sequence is:')
X
p=0:(N-1);
stem(p,abs(X)),grid

The given sequence is x(n): [1 2 2 2 2]

Figure P3: 2 The DFT of the input sequence x(n)

(b) MATLAB Program to Compute FFT of a Periodic or Non-periodic


Signal

54 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

% program 2.3
% to compute FFT using DSP function fft
x=input('The given sequence is x(n): ');
X = fft(x); % calculates the FFT X(k) of the vector x(n)
Xs = fftshift(X); % shifts the vector X(k) in symmetric form
Xsm = abs(Xs) ; % magnitude spectrum
Xsp = angle(Xs) ; % phase spectrum
disp ('The FFT of x(n) is: ');
X;
N=length(X);
P=1:N;
stem(P,X),grid

The given sequence is x(n): [1 2 3 6 5 8]

Figure P3: 3 FFT using DSP function fft


(c) MATLAB Program to Compute FFT

55 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

% program 2.4
% to compute FFT using DSP function fft
n = 0:29;
x1 = cos(2*pi*n/10); % 3 periods
x2 = [x1 x1]; % 6 periods
x3 = [x1 x1 x1]; % 9 periods
N = 2048;
X1 = abs(fft(x1,N));
X2 = abs(fft(x2,N));
X3 = abs(fft(x3,N));
F = (0:N-1)/N;
subplot(3,1,1)
plot(F,X1),title('3 periods');
axis([0 1 0 50]);
subplot(3,1,2)
plot(F,X2);
title('6 periods');
axis([0 1 0 50])
subplot(3,1,3)
plot(F,X3);
title('9 periods')axis([0 1 0 50]);

Figure P3: 4 compute FFT

5. Lab Task:-
1. Try to explain step by step how the programs 1.2, 1.3 and 1.4 work.
2. Modify programs 1.4 to accept any arbitrary user data sequences.
3. Run the programs for different input data sequences.

56 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Lab #4:- Z-Transform Analysis


1. Objectives:
 To analyze signals in the complex frequency domain (z-transform analysis).
 To convert signals in the z-rational form to partial fraction expansion form and vice
versa.
 To plot pole-zero diagrams of the function in complex frequency domain.

2. Background Theory
The z-transform of the signal x(n) is given by

𝑋(𝑧) = ∑ 𝑥(𝑛)𝑧 −𝑛
𝑛=−∞

Where z is the complex variable (z= 𝜎 +jω).


The inverse z-transform can be obtained analytically as
1
𝑥 (𝑛 ) = ∮ 𝑥(𝑧)𝑧 𝑛−1 𝑑𝑧𝑐
2𝜋𝑗
Where C is a counter-clockwise closed path in the z-plane.
(a) Rational Z-transform to partial fraction form:
This technique is usually used, while taking the inverse Z-transform and when the order of H(z)
is high so that it is quite difficult to solve it mathematically.
Change the following transfer function in equation (3.10) to partial fraction form.
5
2 − 6 𝑧 −1
𝐻 (𝑧 ) =
5 1
1 − 6 𝑧 −1 + 6 𝑧 −2

3. Materials Required: Matlab Software and Computers.


(4a) Procedure:
 Use Matlab command residuez.
 Enter numerator coefficients.
 Enter denominator coefficients.
 Produce numerator and denominator coefficients for partial fraction using residuez
command.

57 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

% clc
n=[2 -5/6]; % numerator coefficients
d=[1 -5/6 1/6]; % denominator coefficients
[p,q,r]=residuez(n,d) % produces numerator and denominator

% coefficients for partial fraction expansion are


p=
1.0000
1.0000
q=
0.5000
0.3333
r=
[]
Hence the partial fraction form is
1 1
𝐻 (𝑧 ) = 1 + 1
1− 𝑧 −1 1− 𝑧 −1
2 3

(b) Partial Fraction Form to Rational Form:


This technique is used when it is required to convert partial fraction expression in to its rational
Z-transform. For example, let us bring equation back into its rational form by using MATLAB
command.

(4b) Procedure:
 Enter numerator coefficients i.e. P values.
 Enter denominator coefficients i.e. q values.
 Use residuez command to produce numerator and denominator coefficients for rational
from.

clc
[a,b]=residuez(p,q,r) % produces numerator and denominator
% coefficients for the rational Z-transform

Hence the rational Z-transform is


a=
2.0000 -0.8333

58 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

b=
1.0000 -0.8333 0.1667
Hence the rational Z-transform is
5
2−0.8333𝑧 −1 2− 𝑧 −1
𝐻 (𝑧 ) = = 5
6
1
1−0.8333𝑧 −1+0.1667𝑧 −2 1− 𝑧 −1 + 𝑧 −2
6 6

(c) Pole Zero Diagrams for a Function in Z-domain:


The MATLAB command ‘zplane’ computes and displays the pole-zero diagram of Z-transform.
Let us compute equation (3.10).

clc
n=[2 -5/6]; % numerator coefficients
d=[1 -5/6 1/6]; % denominator coefficients
roots(n) % display zeros
roots (d) % display poles
zplane(n,d) % displays the pole-zero diagram

ans =
0.4167
ans =
0.5000
0.3333

Figure P3: 5 Pole- zero diagram of H (z)

59 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Hence the region of convergence (ROC) of the given Z-transform is |z |> 1/2.

5. Lab Task:-
1. Express the following Z-transform in the partial fraction expansion form, plot its poles and
zeros, and then determine its ROCs.
2𝑧 4 + 16𝑧 3 + 44𝑧 2 + 56𝑧 + 32
𝐺 (𝑧 ) =
3𝑧 4 + 3𝑧 3 − 15𝑧 2 + 18𝑧 − 12
2. Write a MATLAB program to compute and display the poles and zeros, to compute and
display the factored form, and to generate the pole-zero plot of a z-transform that is as ratio of
two polynomials in z−1. Use this program to analyze G(z),
2 + 5𝑧 −1 + 9𝑧 −2 + 5𝑧 −3 + 3𝑧 −4
𝐺 (𝑧 ) =
2 + 45𝑧 −1 + 2𝑧 −2 + 𝑧 −3 + 𝑧 −4

60 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Experiment Part Four: Sampling Process


In this section, the sampling process of continuous-time analog signals in time domain will be
simulated and studied.
 To simulate the sampling process of analog continuous signals in time domain.
 To study the effect of aliasing in time and frequency domain.
Typically, the discrete-time signals are formed periodically by sampling continuous-time signal
xa(t) as
𝑥(𝑛) = 𝑥𝑎(𝑛𝑇𝑠)
Where Ts is the sampling period and fs=1/Ts is the sampling frequency in samples per second.
In periodic sampling, a continuous-time signal is first multiplied by periodic sequences of
impulses,

𝑠𝑎 (𝑡) = ∑ 𝛿(𝑡 − 𝑛𝑇𝑠 )


𝑛=−∞

The sampled signal xs(t) is formed by multiplying sa(t) and xa(t) as


𝑥𝑠 (𝑡) = 𝑥𝑎 (𝑡)𝑠𝑎 (𝑡) = ∑ 𝑥𝑎(𝑛𝑇𝑠) 𝛿(𝑡 − 𝑛𝑇𝑠 )


𝑛=−∞

The Fourier transform of the sampled signal xs(t) is


𝑥𝑠 (𝑗𝛺) = ∑ 𝑥𝑎 (𝑛𝑇𝑠 ) 𝑒 −𝑗𝛺𝑇𝑠


𝑛=−∞

Another expression for xS( jΩ) follows by noting that the Fourier transform of sa(t) is

2𝜋
sa (𝑗𝛺) = ∑ 𝛿(𝛺 − 𝑘𝛺𝑠 )
𝑇𝑠
𝑛=−∞

Where Ωs = 2 /Ts is the sampling frequency.


Therefore,

61 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual


1 1
𝑋𝑠 (𝑗𝛺) = 𝑋𝑎 (𝑗𝛺) ∗ 𝑆𝑎 (𝑗𝛺) = ∑ 𝑋𝑎 (𝛺 − 𝑘𝛺𝑠 )
2𝜋 𝑇𝑠
𝑛=−∞

Finally, the Fourier transform of x(n) is



𝑗𝑤
𝑋(𝑒 ) = ∑ 𝑋𝑎 (𝑛𝑇𝑠 )𝑒 −𝑗𝜔𝑛
𝑛=−∞

𝑗𝑤 )
1 𝑗𝑤 2𝜋𝑘
𝑋(𝑒 = 𝑋𝑠 (𝑗𝛺)|𝛺=𝑤/𝑇𝑠 = ∑ 𝑋𝑎 ( − 𝑗 )
𝑇𝑠 𝑇𝑠 𝑇𝑠
𝑛=−∞

Thus, X(ejω) is a frequency scaled version of ( Ω), with the scaling defined by
𝑤 = 𝛺𝑇𝑠
This scaling, which makes X(ejω) periodic with a period of 2π, is a consequence of time-scaling
that occurs when xs(t) is converted to x(n).
Sampling Theorem: if xa(t) is strictly bandlimited,
𝑋𝑎 (𝑗𝛺)=0 |Ω|>Ω0
Then xa(t) may be uniquely recovered from its samples xa(nT s) if
2𝜋
𝛺= ≥ 2 𝛺0
𝑇𝑠
The frequency Ω is called the Nyquist frequency, and the minimum sampling frequency,
Ω = 2Ω is called the Nyquist rate.
If the signal is sampled below the sampling rate i.e. Ω < 2Ω0 , there is the overlapping of the
shifted replicas of Xa (Ω) . This phenomenon is known as aliasing. When aliasing occurs the
original signal cannot be recovered.
To generate the reconstructed signal ya(t) from x[n], we pass x[n] through an ideal low pass
reconstruction filter that in turn can be implemented according to the bellow equation . h(t) is the
impulse response of the filter.
𝜋𝑡
sin( 𝑇 )
𝑠
ℎ (𝑡 ) = 𝜋𝑡
𝑇𝑠

62 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

If this equation is computed at closely spaced values of t, a plot of ya(t) will resemble a
continuous-time signal.

Lab #1Sampling of Continuous-Time Signal into Discrete-Time Signal


1. Objective:
 To sample continuous-time signal into discrete-time signal.

2. The background theory about sampling process is written at the beginning .

3. Materials Required: Matlab Software and Computers.

4. Procedure:
% p4.1
% Illustration of the Sampling Process
% in the Time Domain
clc;
t = 0:0.0005:1;
f = 13;
xa = cos(2*pi*f*t);
subplot(2,1,1)
plot(t,xa,'LineWidth',1.5);
xlabel('Time, msec');ylabel('Amplitude');title('Continuous-time signal
x_{a}(t)');
axis([0 1 -1.2 1.2])
subplot(2,1,2);
T = 0.1;
n = 0:T:1;
xs = cos(2*pi*f*n);
k = 0:length(n)-1;
stem(k,xs,'r');
xlabel('Time index n');ylabel('Amplitude');title('Discrete-time signal
x[n]');
axis([0 (length(n)-1) -1.2 1.2])

n = 0:T:1;

xs = cos(2*pi*f*n);
k = 0:length(n)-1;

stem(k,xs,'r');

xlabel('Time index n');ylabel('Amplitude');

title('Discrete-time signal x[n]');


axis([0 (length(n)-1) -1.2 1.2])

63 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Figure P4: 1 Sampling of continuous time signal to discrete-time signal.

5.Lab Task:-

1. What function is sampled in program P4.1?

2. Explain the implementation of program P4.1 in detail.

3. Modify and run the program P4.1 for exponential functions.

4. The plots of the continuous-time signal and its sampled version generated by running Program
P4.1 for the following four values of the sampling period are shown below: T = 0.004 msec, T =
0.02 msec, T = 0.15 msec, T = 0.2 msec. What is your observations? Comment your result.

64 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Lab #2 Aliasing Effect in Time Domain


1. Objective:

 To study the aliasing effect of signals during sampling process in time domain.

2. The background theory about sampling process is written at the beginning .

3. Materials Required: Matlab Software and Computers.

4. Procedure:
Here we will generate a continuous-time equivalent ya(t) of the discrete-time signal x[n]
generated in Program P1.1 to investigate the relation between the frequency of the sinusoidal
signal xa(t) and the sampling period. To generate the reconstructed signal ya(t) from x[n], we
pass x[n] through an ideal low pass reconstruction filter that in turn can be implemented
according to equation (4.12) . If this equation is computed at closely spaced values of t, a plot of
ya(t) will resemble a continuous-time signal.

% Program P4.2
% Illustration of Aliasing Effect in the Time Domain
clc;
T = 0.1;
f = 13;
n = (0:T:1)';
xs = cos(2*pi*f*n);
t = linspace(-0.5,1.5,500)';
ya = sinc((1/T)*t(:,ones(size(n))) - (1/T)*n(:,ones(size(t)))')*xs;
plot(n,xs,'bo',t,ya, 'r','Linewidth',1.5);grid;
xlabel('Time, msec');ylabel('Amplitude');
title('Reconstructed continuous-time signal y_{a}(t)');
axis([0 1 -1.2 1.2]);

Figure P4: 2 Aliasing Effect in the Time Domain

65 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

5. Lab Task:-
1) Run Program 4.2 to generate both the discrete-time signal x[n] and its continuous time
equivalent ya(t) and display them.
2) What is the range of t and the value of the time increment in Program 4.2? What is the
range of t in the plot? Change the range of t so as to display the full range ya(t) being
computed in the above program and run Program 4.2 again. Comment on the plot
generated after this change.
3) The plots of the continuous-time signal and its sampled version generated by running
Program P4.2 for the following four values of the sampling period are shown below: T =
0.004 msec, T = 0.02 msec, T = 0.15 msec, T = 0.2 msec. What is your observations?

66 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Lab #3 Effect of Sampling in Frequency Domain: Aliasing Effect in


Frequency
1. Objective:
 To study the effect of aliasing when the signal is represented in frequency domain.

2. The background theory about sampling process is written at the beginning.

3. Materials Required: Matlab Software and Computers.

4. Procedure:
The relation between the continuous-time Fourier transform (CTFT) of an arbitrary
bandlimited Continuous-time signal and the discrete-time Fourier transform (DTFT) of the
discrete-time signal is investigated next. In order to convert a continuous-time signal xa(t) into
an equivalent discrete-time signal x[n], the former must be band-limited in the frequency
domain. To illustrate the effect of sampling in the frequency domain we choose an exponentially
decaying continuous-time signal with a CTFT that is approximately bandlimited.

67 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

% Program 4.3
% Illustration of the Aliasing Effect
% in the Frequency Domain
clf;
t = 0:0.005:10;
xa = 2*t.*exp(-t);
subplot(2,2,1)
plot(t,xa);grid
xlabel('Time, msec');ylabel('Amplitude');
title('Continuous-time signal x_{a}(t)');
subplot(2,2,2)
wa = 0:10/511:10;
ha = freqs(2,[1 2 1],wa);
plot(wa/(2*pi),abs(ha));grid;
xlabel('Frequency, kHz');ylabel('Amplitude');
title('|X_{a}(j\Omega)|');
axis([0 5/pi 0 2]);
subplot(2,2,3)
T = 1;
n = 0:T:10;
xs = 2*n.*exp(-n);
k = 0:length(n)-1;
stem(k,xs);grid;
xlabel('Time index n');ylabel('Amplitude');
title('Discrete-time signal x[n]');
subplot(2,2,4)
wd = 0:pi/255:pi;
hd = freqz(xs,1,wd);
plot(wd/(T*pi), T*abs(hd));grid; xlabel('Frequency,
kHz');ylabel('Amplitude'); title('|X(e^{j\omega})|'); axis([0 1/T 0 2])

Figure P4: 3 Effect of aliasing in frequency domain

68 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

5. Lab Task:-
1. Run Program 4.3 to generate and display both the discrete-time signal and its continuous-time
equivalent, and their respective Fourier transforms. Is there any visible effect of aliasing?
2. Repeat Program 1.3 by increasing the sampling period to 1.5. Is there any visible effect of
aliasing?
2
3. Modify Program 4.3 for the case of 𝑥𝑎 (𝑡) = 𝑒 −𝑖𝜋𝑡 and repeat questions 1 and 2.
4. Time-domain Sampling: Consider the continuous-time signal
𝑒 −2𝑓𝑜 𝑡 , 𝑡 ≥ 0
𝑋𝑎 (𝑡) = {
0, 𝑡 < 0

(a) Compute analytically the spectrum Xa(f) of xa(t).


(b) Compute analytically the spectrum of the signal x(n)=xa(nT), T=1/fs.
(c) Plot the magnitude spectrum | 𝑋𝑎 (𝑓)|for fo=10Hz.
(d) Plot the magnitude spectrum | 𝑋𝑎 (𝑓)|for fs=10, 20, 40, and 100Hz.
(e) Explain the results obtained in part (d) in terms of the aliasing effect.

69 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Lab #4 DECIMATION & INTERPOLATION


1. Objective
 To investigate using MATLAB the operations of the up-sampler and the down-sampler
both in the time domain and in the frequency domain.
2. Background theory
The digital signal processing structures discussed so far belong to the class of single-rate
systems as the sampling rates at the input and the output and all internal nodes are the same. There
are applications where it is necessary and often convenient to have unequal rates of sampling at
various parts of the system including the input and the output. In this laboratory exercise you will
investigate first using MATLAB the properties of the up-sampler and the down-sampler, the two
basic components of a multi-rate system. You will then investigate their use in designing more
complex systems, such as interpolators and decimators, and filter banks.
3. Materials Required: Matlab Software and Computers.
4. Procedure:

a) Input-Output Relations in the Time-Domain:-


Program 4.4
% Illustration of the Sampling Process
% in the Time Domain
clf;
t = 0:0.0005:1;
f = 13;
xa = cos(2*pi*f*t);
subplot(2,1,1)
plot(t,xa);grid
xlabel('Time, msec');ylabel('Amplitude');
title('Continuous-time signal x_{a}(t)');
axis([0 1 -1.2 1.2])
subplot(2,1,2);
T = 0.1;
n = 0:T:1;
xs = cos(2*pi*f*n);
k = 0:length(n)-1;
stem(k,xs); grid
xlabel('Time index n’);ylabel(’Amplitude');
title('Discrete-time signal x[n]');
axis([0 (length(n)-1) -1.2 1.2])

70 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Figure P4: 4 Illustration of the Sampling Process

%ProgramP4_5
%Illustration of Down-Sampling by an Integer Factor
%
close all; clear all; clc
n = 0:49;
m = 0:50*3-1;
x = sin(2*pi*0.042*m);
y = x([1:3:length(x)]);
subplot(2,1,1)
stem(n,x(1:50));axis([0 50 -1.2 1.2]);
title('InputSequence');
xlabel('Timeindexn');
ylabel('Amplitude');
subplot(2,1,2)
stem(n,y);axis([0 50 -1.2 1.2]);
title('OutputSequence');xlabel('Timeindexn');ylabel('Amplitude');

Figure P4: 5 Illustration of Down-Sampling by an Integer Factor


The Signal Processing Toolbox includes three M-functions which can be employed to design and
implement an interpolator or a decimator. The three M-functions are decimate, interp, and

71 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

resample. Each function is available with several options. In this section you will study the
decimation and interpolation operation using these functions.
b) Decimator Design and Implementation
Program P10_3 illustrates the use of the M-function decimate in the design and implementation
of a decimator with an integer-valued decimation factor M. In the option utilized in this program,
decimate designs and uses a low pass decimation filter with a stopband edge.
%ProgramP4_6
%Illustration of Decimation Process
clear all; close all; clc
M=input('Down-samplingfactor=');
n=0:99;
x=sin(2*pi*0.043*n)+sin(2*pi*0.031*n);
y=decimate(x,M,'fir');
subplot(2,1,1);
stem(n,x(1:100));
title('InputSequence');
xlabel('Timeindexn');ylabel('Amplitude');
subplot(2,1,2);
m=0:(100/M)-1;
stem(m,y(1:100/M));
title('OutputSequence');xlabel('Timeindexn');ylabel('Amplitude');

Down-samplingfactor=5

Figure P4: 6 Illustration of Decimation Process

c) Interpolator Design and Implementation


Program P4_7 illustrates the use of the M-function interp in the design and implementation
of an interpolator with an integer-valued interpolation factor L. interp designs and uses a
lowpass interpolation filter with a stopband edge satisfying.
𝐿, |𝜔| ≤ 𝜔𝑐 /𝐿
|𝐻(𝑒 𝑗𝜔 )| = {
0, 𝜋/𝐿 ≤ |𝜔| ≤ 𝜋

72 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

%ProgramP4-7
%Illustration of Interpolation Process
%
clear all; close all; clc
L=input('Up-samplingfactor=');
% Generate the input sequence
n=0:49;
x=sin(2*pi*0.043*n)+sin(2*pi*0.031*n);
% Generate the interpolated output sequence
y=interp(x,L);
% Plot the input and the output sequences
subplot(2,1,1);
stem(n,x(1:50));
title('InputSequence');xlabel('Timeindexn');ylabel('Amplitude');
subplot(2,1,2);
m=0:(50*L)-1;
stem(m,y(1:50*L));
title('OutputSequence');xlabel('Timeindexn');ylabel('Amplitude');

Up-samplingfactor=2

Figure P4: 7 Illustration of Interpolation Process


d) Fractional-Rate Sampling Rate Alteration
Program P4_8 illustrates the use of the M-function resample in the design and implementation
of an interpolator with a fractional-rate interpolation factor L/M. Resample designs and uses a
lowpass interpolation filter with a stopband edge.

73 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

% Program P4-8
% Illustration of Sampling Rate Alteration by
% a Ratio of Two Integers
%
close all; clear all; clc
L=input('Up-samplingfactor=');
M=input('Down-samplingfactor=');
n=0:29;
x=sin(2*pi*0.43*n)+sin(2*pi*0.31*n);
y=resample(x,L,M);
subplot(2,1,1);
stem(n,x(1:30));axis([0 29 -2.2 2.2]);
title('InputSequence');
xlabel('Timeindexn');ylabel('Amplitude');
subplot(2,1,2);
m=0:(30*L/M)-1;
stem(m,y(1:30*L/M));axis([0 (30*L/M)-1 -2.2 2.2]);
title('OutputSequence');
xlabel('Timeindexn');ylabel('Amplitude');

Up-samplingfactor=1
Down-samplingfactor=2

Figure P4: 8 Illustration of Sampling Rate Alteration


5. Lab Task
1. What is the angular frequency in radians of the sinusoidal sequence in Program P4_1?
What is its length? What is the up-sampling factor L?
2. How is the up-sampling operation implemented in Program P4_1?
3. Modify Program P4_1 to study the operation of an up-sampler on a ramp sequence.
4. Program P4_2 can be used to study the operation of a down-sampler
5. What is the angular frequency in radians of the sinusoidal sequence Program P4_2?
What is its length? What is the down-sampling factor M?

74 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

6. How is the down-sampling operation implemented in Program P4_2?


7. What are the frequencies of the two sinusoidal sequences forming the input sequence in
Program P4_3? What is the length of the input?
8. What are the type and order of the decimation filter?
9. Run Program P4_3 for M = 4 and comment on the results.
10. Change the frequencies of the two sinusoidal sequences in Program P4_3 in the input
signal to 0.045 and 0.029, and the length of the input to 120. Run the modified Program
P10 5 for M = 3. Comment on your results.

75 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Experiment Part Five: Digital Filter Design


In this part of the experiment, the processes of digital filter designs using Matlab will be
simulated, and studied. The two basic digital filter design methods to be characterized and
studied are the following: Infinite impulse response (IIR) based on impulse invariance method
and bilinear transformation method, and design of finite impulse response (FIR) filters by
windowing methods.
It is clear that digital filters are the fundamental building blocks in digital signal processing since
they play a paramount role in removing and reducing unwanted noise from the desired signal.
Lab #1 Infinite Impulse Response (IIR) Filter Design Based on the
Impulse Invariance Method and Bilinear Transformation
1. Objectives:
 To study and simulate design process of IIR digital filters based on impulse invariance
method.
 To study and simulate digital filter design process of IIR design based on bilinear
transformation.
2. Background Theory
(a) Infinite Impulse Response (IIR) Filter Design Based on the Impulse Invariance
Method
For transforming an analogue continuous filter to an IIR digital filter using the Impulse
Invariance method, the digital IIR filter H(z) is obtained by multiplying the sampling period Ts
by the z-transform of H(s).
Suppose the transfer function or system function of analog continuous filter is H(s). H(s) is
brained as follows:
1. Transform the specifications of digital filter into continuous analog filter
specifications. 2. Obtain the order of the filter using well known analog filter types such
as Butterworth, Chebyshev, elliptic and Linear Phase filters from the transformed
specifications.
3. Obtain the poles Sk of the analog filter using of Butterworth, Chebyshev, elliptic or
Bessel types from the standard formula table.
4. Determine the transfer function H(s) of analog filter as follows from its poles:
𝑁
−𝑠𝑘
𝐻 (𝑠 ) = ∏
𝑠 − 𝑠𝑘
𝑘=1

5. Express H(s) in partial fraction as follows:


𝑁
𝐴𝐾
𝐻 (𝑠 ) = ∑
𝑠 − 𝑠𝑘
𝑘=1

76 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Once you obtain partial fraction coefficients Ak, determine the transfer function of IIR digital
filter H(z) as follows
𝑁
𝑇𝑠 𝐴𝐾
𝐻 (𝑧 ) = ∑
1 − 𝑒 𝑠𝑘𝑇𝑠 𝑧 −1
𝑘=1

Ts is the sampling time and it should be chosen appropriately.


MATLAB contains the function impinvar that is used to produce IIR digital filters based on the
impulse invariance transformation.
The syntax is
[numz,denz] = impinvar(nums,dens,fs)
This function creates a digital filter with numerator and denominator coefficients numz and denz
respectively, whose impulse response is equal to the impulse response of the analog filter with
coefficients nums and dens and a sampling frequency fs. The nums and dens coefficients will be
scaled by fs. If you don't specify fs, it defaults to 1 Hz.
The syntax
[numz,denz] = impinvar(nums,dens,fs,tol) uses the tolerance tol for grouping repeated poles
together. Default value is 0.001, i.e., 0.1%. The repeated pole case works, but is limited by the
ability of the function roots to factor such polynomial.
(b) Infinite Impulse Response (IIR) Filter Design Based on Bilinear Transformation Method
For transforming an analogue continuous filter to an IIR digital filter using the Bilinear method,
the digital IIR filter transfer function H(z) is obtained as
𝐻 (𝑧) = 𝐻(𝑠)| 2 1+𝑧 −1
𝑠= ( )
𝑇𝑠 1+𝑧 −1

2 1+𝑧 −1
In other words, evaluate H(s) at 𝑠 = 𝑇 (1+𝑧 −1 )
𝑠

Where H(s) is the transfer function of analog filter.


The frequency transformations are the following:
2
𝛺= 𝑡𝑎𝑛(𝑤/2)
𝑇𝑠
Where Ω is analog filter’s frequency and ω is the frequency of digital filter.

77 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

With higher order filter the calculations get very complex. MATLAB contains the function
bilinear and is used to obtain the transfer function H(z) that approximates the analogue transfer
function H(s).
Within the syntax
[zz,pz,kz] = bilinear(zs,ps,ks,fs)
the function converts the s-domain transfer function specified by zs, ps and ks to a z-transform
discrete equivalent obtained from the Bilinear transformation where column vectors zs and ps
specify the zeros and poles, scalar ks specifies the gain and fs is the sampling frequency in Hz.
With the syntax
[numz,denz] = bilinear(nums,dens,fs)
where nums and dens are row vectors containing the numerator and denominator coefficients of
the analogue transfer function H(s), in descending powers of s, the function transforms H(s) to
the z-transform coefficients numz and denz of the IIR digital filter transfer function H(z). The
syntax [Az,Bz,Cz,Dz] = bilinear(As,Bs,Cs,Ds,fs) is a state-space version of the bilinear
transform. Each syntax for the bilinear function accepts an optional additional input argument
that specifies prewarping.
The syntax
[zz,pz,kz] = bilinear (zs,ps,ks,fs,fprewarp)
Applies prewarping before the bilinear transformation so that the frequency responses before and
after mapping match exactly at frequency point fprewarp. The matching point fprewarp is
specified in Hz.
(c) IIR Filter Design Using MATLAB
(i) From the Analogue Prototype to the IIR Digital Filter
In this way of design, analog filter of any type can be used and transformed into IIR digital
equivalent filters using impulse invariance and bilinear transformations. The two transform
MATLAB functions are impinvar and bilinear.
Table P5: 0-1 Analogue Transformation Functions
Frequency Transformation Transformation Function
Lowpass to lowpass [nums,dens]=lp2lp(nums,dens,W0)
[As,Bs,Ds]=lp2lp(As,Bs,Cs,Ds,W0)
Lowpass to highpass [nums,dens]=lp2hp(nums,dens,W0)
[As,Bs,Ds]=lp2hp(As,Bs,Cs,Ds,W0)
Lowpass to bandpass [nums,dens]=lp2bp(nums,dens,W0,BW)
[As,Bs,Ds]=lp2bp(As,Bs,Cs,Ds,W0,BW)

78 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Lowpass to bandstop [nums,dens]=lp2bs(nums,dens,W0, BW)


[As,Bs,Ds]=lp2bs(As,Bs,Cs,Ds,W0,BW)
Table P5: 0-2 Transformation from Analogue to Digital Filters

Frequency Transformation Transformation Function


Bilinear [zz,pz,kz]=bilinear(zs,ps,ks,fs)
[zz,pz,kz]=bilinear(zs,ps,ks,fs,fprewarp)
Numz,denz]=bilinear(nums,dens,fs,fprewarp)
Impulse Invariance [numz,denz]=impinvar(nums,dens,fs)
[numz,denz]=impinvar(nums,dens,fs,tol)

(ii) Direct Design


MATLAB has many functions that will help us design digital IIR filters directly. These
MATLAB functions are listed in Table below. In Table below, n is the order of the desired IIR
digital filter. wn is the single-valued cut-off frequency of the lowpass filter. In case of a bandpass
or a bandstop digital filter, wn is a vector of two values: the lower and the upper cut-
offfrequencies. Rp and Rs are the passband and the stopband attenuation ripples in decibels.
ftype is the type of the desired filter. We use stop for stopband digital filters and high for
highpass digital filters. If you do not specify the ftype the default is lowpass. The wn
frequency/frequencies must be normalized to the range [0 1]. The range [0 1] corresponds to the
range [0 π]. zz, pz and kz are the zeros, poles and gain constants.
Table P5: 0-3 Direct Design Matlab Functions for IIR Filters

Filter Type Analogue Prototype Function


Butterworth [numz,denz]=butter(n,wn,’ftype’)
[zz,pz,kz]=butter(n,wn,’ftype’)
Chebyshev Type I [numz,denz]=cheby1(n,Rp,wn,’ftype’)
[zz,pz,kz]=cbeby1(n,Rp,wn,’ftype’)
Chebyshev Type II [numz,denz]=cheby2(n,Rp,wn,’ftype’)
[zz,pz,kz]=cbeby2(n,Rp,wn,’ftype’)
Elliptic [numz,denz]=ellip(n,Rp,wn,’ftype’)
[zz,pz,kz]=ellip(n,Rp,wn,’ftype’)

3. Materials Required: Matlab Software and Computers.

4. Procedure:
Start with a prototype second-order analogue lowpass Butterworth filter.
 Convert it to a second-order lowpass analogue Butterworth filter with cut-off frequency
of 200π rad/sec.

79 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

 Use the Bilinear and the Impulse Invariance transformations to get H(z), the transfer
function of the IIR digital filter.
 IIR Impulse invariance and bilinear transformation using direct method:

IIR Impulse invariance and bilinear transformation

a. using direct method:


% program 5.1
%we will generate the lowpass analogue Butterworth prototypes first
for n=2:4:6% N=2 and 6
[zs, ps, ks]=buttap(n);
% put in terms of transfer function
[nums, dens]=zp2tf(zs,ps,ks);
%next we transform to analogue lowpass with cut-off 200pi
[nums, dens]=lp2lp(nums, dens, 200*pi);
%Then we use the Bilinear transform first
fs=1000;
[numzb, denzb]=bilinear(nums, dens, fs);
%then we use the Impulse Invariance method
[numzi, denzi]=impinvar(nums, dens, fs);
%Now we plot the magnitude response of the IIR digital filter
%using the two transformation methods
[Hb, fb]=freqz(numzb, denzb);
[Hi, fi]=freqz(numzi, denzi);
subplot(2,1,1);
plot(fb/pi, abs(Hb));
hold on;
axis([0 1 0 1]);
title('The Magnitude response using the Bilinear transform');
subplot(2,1,2);
plot(fi/pi, abs(Hi));
hold on;
end
title('The Magnitude response using the Invariance transform');
xlabel('frequency in pi units');
axis([0 1 0 1]);

80 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Figure P5. 1 Magnitude responses of IIR filter using impulse invariance and bilinear
transformation by method I.

Using the Bilinear transformation method, the numerator and denominator coefficients of the
transfer function of the IIR digital second-order filter are
numzb = 0.0640 0.1279 0.0640
denzb = 1.0000 –1.1683 0.4241 and the corresponding IIR filter transfer function is
0.0640𝑧 2 + 0.127𝑧 + 0.0640
𝐻 (𝑧 ) =
𝑧 2 − 1.1683𝑧 + 0.4241
Using the Impulse Invariance method, the numerator and denominator coefficients of the transfer
function of the IIR digital second-order filter are
numzi = 0 0.2449
denzi = 1.0000 –1.1580 0.4112 and the corresponding transfer function of the IIR filter is
0.2449
𝐻 (𝑧 ) =
𝑧 2 − 1.1580𝑧 + 0.4112
b. Method II: Direct Design
 Use the MATLAB function : [numz, denz] = butter (n,wn,’ftype’) without ftype, since
the default type is lowpass, and with wn = 0.2 and n= 2 and 6.
 Use wn = 0.2 as the normalized frequency because the equivalent cut-off digital
frequency is 200 π(1/1000) = 0.2π.
 Normalize, using MATLAB, the digital frequencies to belong to the digital frequency
interval [0 1].

81 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

The following MATLAB script will accomplish the design:

% program 5.2
clf
for n=2:4:6
[numz, denz]=butter(n, 0.2)%the normalized digital frequency is 0.2
%Now we plot the magnitude response of the IIR digital filter
%using the two transformation methods
[H, f]=freqz(numz, denz);
plot(f/pi, abs(H));hold on;
end
title('The Magnitude response using the function butter');
xlabel('frequency in pi units');
axis([0 1 0 1]);

In this case, coefficients of the numerator and the denominator of the transfer function for N = 2
are
numz = 0.0003 0.0020 0.0051 0.0068 0.0051 0.0020 0.0003
denz =1.0000 -3.5794 5.6587 -4.9654 2.5295 -0.7053 0.0838

Figure P5. 2 The magnitude responses of IIR filter using impulse invariance and bilinear
transformation by direct method
The corresponding digital IIR transfer functions are
0.0657𝑧 2 +0.1349𝑧+006557
𝐻 (𝑧 ) = 𝑧 2−1.143𝑧+0.4128

5. Lab Task:-
1. Consider the input continuous signal x(t)=2+ cos(10t) Design a digital second- and sixth-
order IIR low pass Butterworth filter to attenuate the cos(10t) term and pass the dc 2 term.
2. We are trying to eliminate all frequencies that are lower than 100 Hz from the incoming
signal x(t). If the sampling frequency is 1000 Hz, design a digital high pass IIR filter to

82 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

accomplish this task. Allow 1 dB for the maximum passband attenuation and 40 dB for the
minimum stopband attenuation.
3. We are interested in passing the term sin(1500t) from the input signal x(t) = sin(10t)
+sin(1500t) + sin(5000t). Design a bandpass IIR Cheby1 digital filter to accomplish this
task. Plot the input signal along with the magnitude response of the filter and its output.
4. We are trying to suppress all frequencies that are higher than 200 Hz and below 400 Hz from
the incoming signal x(t). If the sampling frequency is 2000 Hz (the minimum sampling
should be 800 Hz but sampling at higher rates can give better results), design a digital
bandstop IIR filter to accomplish this task. Use the Cheby2 filter in the design. Use 40 dB for
the minimum allowable ripple in the stopband.

Note: In all questions you have to use Matlab to design the required IIR filters

83 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Lab #2 Filter Design by Matlab tools


1. objective
 To simulate and study the design process of FIR digital filters by using windowing
methods.

2. Background Theory
The general input/output equation of a causal FIR filter in the Z-domain is given by:

𝑌(𝑧) = 𝐻 (𝑧)𝑋(𝑧) = (ℎ(0) + ℎ(1)𝑧 −1 + ⋯ + ℎ(𝑁 − 1)𝑧 (𝑁−1) )𝑋(𝑧)


Here the time origin is chosen equal to zero. In the time domain the equation becomes a
convolution:
𝑦( 𝑛 ) = ℎ ( 0 ) 𝑥 ( 𝑛 ) + ℎ ( 1 ) 𝑥 ( 𝑛 − 1 ) + ⋯ + ℎ ( 𝑁 − 1 ) 𝑥 ( 𝑁 − 1 )
However, because in Matlab the index of a vector needs always starts from sample 1 (not 0), the
first filter coefficient h(0) will be at position 1 its MatLab vector representation hcoef, causing
the MatLab representation of the filtering equation to become:
𝑦(𝑛) = ℎ𝑐𝑜𝑒𝑓(0)𝑥(0) + ℎ𝑐𝑜𝑒𝑓(1)𝑥(𝑛 − 1) … + ℎ𝑐𝑜𝑒𝑓(𝑁 − 1)𝑥(𝑁 − 1)

3. Materials Required: Matlab Software and Computers.

4. Procedure:
Matlab includes several graphical tools to generate filters. The Signal Processing Tool
(sptool) is an interactive GUI for digital signal processing used to not only analyse signals but
also to design and to analyse filters. You can accomplish these tasks using four GUIs that you
access from within sptool:
 The Signal Browser is for analyzing signals. You can also play signals using your
computer’s audio hardware.
 The Filter Design and Analysis Tool (fdatool) is available for designing or editing FIR
and IIR digital filters. Most filter design methods available at the command line are also
available in fdatool. Additionally, you can use fdatool to design a filter by using the
Pole/Zero Editor to graphically place poles and zeros on the z-plane.
 The Filter Visualization Tool (fvtool) is the proper tool to analyses filter, showing its
characteristics.
 The Spectrum Viewer is for spectral analysis. You can use Signal Processing Toolbox
spectral estimation methods to estimate the power spectral density of a signal.

84 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Design a digital FIR lowpass filter by using matlab’s tools few steps.
1) Type sptool at the MATLAB command prompt:

Figure P5. 3 SPTool startup

2) select ’New Design’.

Figure P5. 4 Filter Designer

85 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

The GUI has three main regions:


 The Current Filter Information region
 The Filter Display region and
 The Design panel
The upper half of the GUI displays information on filter specifications and responses for the
current filter. The Current Filter Information region, in the upper left, displays filter properties,
namely the filter structure, order, number of sections used and whether the filter is stable or not.
It also provides access to the Filter manager for working with multiple filters.
The Filter Display region, in the upper right, displays various filter responses, such as,
magnitude response, group delay and filter coefficients.
The lower half of the GUI is the interactive portion of Filter Designer. The Design Panel, in
the lower half is where you define your filter specifications. It controls what is displayed in the
other two upper regions. Other panels can be displayed in the lower half by using the sidebar
buttons.
3) Designing a Filter
Design a digital FIR lowpass filter using the mentioned Matlab’s tool called sptool with the
following specifications:
 Passband cutoff frequency: fp = 2 kHz
 Stopband cutoff frequency: fs = 3 kHz
 Passband Ripple: Rp = 0.25 dB
 Stopband attenuation: Rs = 50 dB
 Sampling frequency: fs = 20 kHz

Figure P5. 5 Designing filter

86 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

4) When you finish, click on the ’Design Filter’ button.

Figure P5. 6 show low pass filter

5) To access the filter coefficients of the designed filter go to the sptool window and select
’Export’ from the ’File Menu’.

Figure P5. 7 Export way


6) Next highlight the designed filter and use the ’Export to Workspace’ button to make the
filter parameters accessible on the workspace. The filter parameters are stored under the
variable name Num in an object-oriented manner. Check the Matlab Workspace.
7) Filter response
In command windows:- Type >> stem(Num)

87 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

Figure P5. 8 show a low pass filter

5. Lab Task:-
1) Design a digital FIR bandpass filter with the following specifications:
Passband: 8-12 kHz
Stopband Ripple: Rs = 80db
Passband Ripple: Rp = 0.001
Transition width: 3kHz
Sampling frequency: fs = 44.1 kHz
Obtain the filter coefficients and frequency response for the above FIR using the Blackman
window method.
2. Design a digital FIR filter for telephony, Specifications for the Transmit Filter:
• attennuation band 1: 0-150 Hz, min. 30 dB attenuation
• pass band: 340 - 3400 Hz, riple max. 2 dB
• attennuation band 2: above 4600 Hz, min. 40 dB attenuation

88 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN


Digital Signal Processing Lab Manual

89 | P a g e DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING by SEREKEBIRHAN

You might also like