Daffodil International University

Department of Electrical and Electronic Engineering

LABORATORY MANUAL FOR

Digital Signal Processing Laboratory

Course No: EEE 322

Course Title: Digital Signal Processing Laboratory
For the students of
Department of Electrical and Electronic Engineering
 Daffodil International University
Dept. of Electrical and Electronic Engineering

EEE 322: Digital Signal Processing Laboratory

1-A. Introduction to MATLAB 1 – 13

1-B. Script Files in MATLAB 14 – 26

2. Implementation of Discrete-time Signals 27 – 34

3. Sampling & Reconstruction of Analog Signals 35 – 40

4. Quantization of Discrete-time Signals 41 – 46

5. Time-Domain Analysis 47 – 53

6. z-Transform and its Application 54 – 61

7. Frequency-Domain Analysis 62 – 70

Revised and updated by –

Saugato Rahman Dhruba
Lecturer, Department of EEE
Daffodil International University
 Daffodil International University
Department of Electrical and Electronic Engineering
EEE 322: Digital Signal Processing Laboratory

EXPERIMENT NO: 01-A

NAME OF THE EXPERIMENT: INTRODUCTION TO MATLAB

The easiest and best way to learn MATLAB is to use MATLAB

What is MATLAB
MATLAB stands for MATrix LABoratory. MATLAB is a software package for computation in engineering, science
and applied mathematics. It offers a powerful programming language, excellent graphics and a wide range of expert
knowledge. MATLAB is published by and a trademark of The MathWorks Inc.

Standard Windows
When the MATLAB program starts, the following window opens which contains four smaller windows. These
windows are default and can be switched on or off if necessary.

Command Current Directory & Command

History window Workspace window window

Figure 1.1: The default view of MATLAB desktop

1
 Four Most Used MATLAB Windows
Window Name What Happens in the Window
Command Window This is where you type and execute commands
This shows current variables and allows to edit variables by opening array
Workspace Window
editor (double click), to load variables from files, and to clear variables
This shows current directory and MATLAB files in current folder, provides
Current Directory Window
with a handy way to change folders and to load files
This shows previously executed commands.
History Window
Commands can be re-executed by double-clicking on them

Working on the Command Window

(A) Using MATLAB as a Calculator
You can use the command window for various calculations as like a calculator: Type the expression containing
numbers, parenthesis and mathematical operations. Press the Enter button to get the result.
Example 1
>> 2^sqrt(9)
ans =
8
Example 2
>> a = 1.2;
>> b = 2.3;
>> c = 4.5;
>> d = 4;
>> a^3 + sqrt(b*d) - 4*c
ans =
-13.23887
Note the semicolon (;) after each variable assignment. If you omit the semicolon, then MATLAB echoes back
on the screen the variable value.

(B) Assigning Values to Variables

To create a variable, just use it on the left hand side of an equal sign (=). The following example shows how to
assign values to variables x, y and z.
Example 3
>> x = 5
x =
5
>> y = log(2)
y =

2
 0.6931
>> z = 10;
>>
In the last case, nothing is displayed. Why? The answer is – because of the use of semicolon (;). Here, a value
of 10 is stored in the variable z, but it is not displayed.

Note:
(1) You can use the „who‟ command to list the currently active variables. For example 3, this results:
>>who
Your variables are:
xyz
(2) You can use the „clear‟ command to delete variables from computer memory:
>>clear x
>>x
Undefined function or variable 'x'.
(3) You can suppress intermediate outputs with semicolon:
MATLAB will print the result of every assignment operation unless the expression on the right hand side is
terminated with a semicolon (See Example 3).

(C) Useful Commands for managing variables

Command Outcome
clear Removes all variables from the memory
clear x y z Removes only variables x, y, z from the memory
who Displays a list the variables currently in the memory
Displays a list the variables currently in the memory and their size, together
whos
with information about their bytes and class

(D) Commands to ensure fairness and readability

Command Outcome
Clears the contents on the command window, resulting a blank working
clc
environment
If used before any line or statement, program ignores the line or statement.
%
Used to insert any comment in the program
Clf Clears the contents on the figure windows and prepare it for a new plot

3
(E) Predefined Variables
A number of frequently used variables are already defined when MATLAB is started. Some of the predefined
variables are:

Predefined Variables Meaning

Value of the last expression that has not been assigned to a specific variable.
ans If the user does not assign the value of an expression to a variable, MATLAB
automatically stores the result in ans
pi Value of the number
eps Smallest difference between two numbers , -
inf Used for infinity ( )
i or j Defined as √ . Used in complex numbers (or variables) operations
Stands for Not-a-Number. Used when MATLAB cannot determine a valid
NaN
numeric value. Example: 0/0

(F) Complex Numbers

MATLAB also supports complex numbers and complex operations. The imaginary number is denoted with the
symbol i or j, assuming that you did not use these symbols anywhere else in your program (very important!).
Example 4
>> z = 3 + 4i
z =
3.0000 + 4.0000i
Note that, you do not need the „*‟ after 4. Here are some useful functions for complex operations:
Example
Function Description
(z is from example 4)
>> conj(z)
conj(z) Computes the complex conjugate of z ans =
3.0000 - 4.0000i
>> real(z)
real(z) Computes the real part of z ans =
3
>> imag(z)
imag(z) Computes the imaginary part of z ans =
4
>> abs(z)
abs(z) Computes the magnitude of z ans =
5
>> angle(z)
angle(z) Computes the phase of z (in radians) ans =
0.9273

4
(G) Elementary Math functions used in MATLAB

Function Description Example

>> sqrt(16)
sqrt(x) Square root of x, i.e., √ ans =
4
>> exp(3)
exp(x) Exponential value of x, i.e., ans =
20.0855
>> abs(-4)
ans =
Absolute value of x, i.e., | | 4
abs(x)
(Both real and complex numbers) >> abs(3 + 4j)
ans =
5
>> log(1000)
log(x) Natural logarithm (Base e) of x, i.e., ( ) ans =
6.9078
>> log10(1000)
log10(x) 10-base logarithm of x, i.e., ( ) ans =
3
>> factorial(5)
The factorial function of
factorial(x) ans =
a positive integer x, i.e.,
120
>> sin(pi/6)
sin(x) Sine of angle x (x in radians) ans =
0.5000
>> cos(pi/6)
cos(x) Cosine of angle x (x in radians) ans =
0.866
>> tan(pi/6)
tan(x) Tangent of angle x (x in radians) ans =
0.5774
>> cot(pi/6)
cot(x) Cotangent of angle x (x in radians) ans =
1.7321
>> asin(sqrt(3)/2)
asin(x) Inverse sine of x (in radians) ans =
1.0472
>> acos(sqrt(3)/2)
acos(x) Inverse cosine of x (in radians) ans =
0.5236
>> atan(1/2)
atan(x) Inverse tangent of x (in radians) ans =
0.4636
>> round(17/5)
ans =
3
round(x) Round the value of x to the nearest integer
>> round(18/5)
ans =
4
>> fix(13/5)
fix(x) Round the value of x towards zero ans =
2

5
 >> ceil(11/5)
Round the value of x towards
ceil(x) ans =
positive infinity ( )
3
>> floor(9/4)
Round the value of x towards
floor(x) ans =
negative infinity ( )
2
>> rem(13,5)
Returns the remainder after
rem(x,y) ans =
x is divided by y
3
>> mod(13,5)
mod(x, y) Returnsthe modulus after division of x by y ans =
3
Sign of a variable x, i.e., >> sign(-5)
sign(x) ans =
( ) 2 -1

(H) Creating Vectors and Matrices

The array or matrix is a fundamental format used to store and manipulate data. Unlike C programming
language, array index starts from 1 in MATLAB instead of 0. To create a matrix in MATLAB, use the square
brackets, [ ]. Numbers or variables inside the brackets can be separated by commas, spaces or semicolons,
where commas and spaces are row separators and semicolons are column separators. Try the following
example.
Example 5
>> a = [2, 5, 7]
a =
2 5 7
>> b = [3, 4, a] % Example of using a variable
b =
3 4 2 5 7
>> c = [2 3;4 5]
c =
2 3
4 5

(I) Colon Operator (:)

MATLAB colon operator is a compact way to create vectors.
Example 6
>> A = 1:0.1:1.5
A =
1.0000 1.1000 1.2000 1.3000 1.4000 1.5000
Where, 0.1 is the increment.
If the increment is omitted, it is assumed to be 1.

6
 Example 7
>> B = 1:5
B =
1 2 3 4 5

(J) MATLAB built-in Matrix Functions

Commands Description Example
>> zeros(2, 3)
ans =
zeros(m, n) A matrix of dimension m × n having all elements 0 0 0 0
0 0 0
>> ones(2, 3)
ans =
ones(m, n) A matrix of dimension m × n having all elements 1 1 1 1
1 1 1
>> zeros(2)
ans =
zeros(n) A matrix of dimension n × n having all elements 0 0 0
0 0
>> ones(2)
ans =
ones(n) A matrix of dimension n × n having all elements 1 1 1
1 1
>> eye(3)
ans =
eye(n) An identity matrix of dimension n×n 1 0 0
0 1 0
0 0 1

(K) Determining the size of a matrix

To determine the size of a matrix, the „size(x)‟ function is used.
Example 7
>> x = [1 2 3; 4 5 6; 7 8 9]
x =
1 2 3
4 5 6
7 8 9
>> size(x)
ans =
3 3
Here, the first element is the number of rows and the second one is the number of columns.

7
(L) Addressing Vectors and Matrices
Here are some commands for addressing matrices or vectors. For illustration, let –

[ ]
In MATLAB –
>> A = [1 3 5 7 9 11;…
2 4 6 8 10 12;…
3 6 9 12 15 18;…
4 8 12 16 20 24;…
5 10 15 20 25 30];

Command Meaning Example

>> B = A(:, 3)
B =
5
Refers to the elements in all the rows of
A(:, n) 6
column n of the matrix A
9
12
15
>> C = A(2, :)
Refers to the elements in all the columns
A(n, :) C =
of row n of the matrix A
2 4 6 8 10 12
>> D = A(:, 2:4)
D =
Refers to the elements in all the 3 5 7
A(:, m:n) rowsbetween columns m and n of the 4 6 8
matrix A 6 9 12
8 12 16
10 15 20
>> E = A(2:4, :)
Refers to the elements in all the E =
A(m:n, :) columnsbetween rows m and n of the 2 4 6 8 10 12
matrix A 3 6 9 12 15 18
4 8 12 16 20 24
>> F = A(1:3, 2:4)
F =
Refers to the elements in rows m through n
A(m:n, p:q) 3 5 7
and columns p through q of the matrix A
4 6 8
6 9 12

(M) Generation of Random Numbers

Command Description Example

>> rand
Generates a single random number
rand ans =
between 0 and 1
0.2311

8
 >> a = rand(1,3)
Generates an n elements row vector of
rand(1, n) a =
random numbers between 0 and 1
0.6068 0.4860 0.8913
>> b = rand(2)
Generates an nn matrix of b =
rand(n)
random numbers between 0 and 1 0.4568 0.4447
0.0158 0.6154
>> c = rand(3,2)
c =
Generates an mn matrix of 0.8147 0.9134
rand(m, n)
random numbers between 0 and 1 0.9058 0.6324
0.1270 0.0975
Generates a row vector with n elements >> randperm(5)
randperm(n) that are random permutation ans =
of integers 1 through n 5 4 2 1 3

(N) Matrix Operations

(1) Matrix Addition and Subtraction
Let, two matrices A and B are defined as –

[ ] [ ]
Then, the sum of A and B is –
[ ]
And, the difference of A and B is –
[ ]

Example: Let,
0 1 0 1
Then, the sum –
0 1 0 1
And, the difference –
0 1 0 1

Example 8
>> A = [1 2 3;4 5 6];
>> B = [8 9 10;11 12 13];
>> A + B
ans =
9 11 13
15 17 19
>> A - B
ans =
-7 -7 -7
-7 -7 -7

9
(2) Matrix Multiplication
Let, two matrices A and B are defined as –

[ ] [ ]

Be sure that, the number of columns of Matrix A is equal to the number of rows of Matrix B. Otherwise,
multiplication of A and B is not defined. The multiplication of A and B is then –

[ ][ ] [ ]

Example: Let,

[ ] [ ]

Then –

[ ][ ] [ ] [ ]

Example 9
>> A = [1 4 3; 2 6 1; 5 2 8];
>> B = [5 4; 1 3; 2 6];
>> A*B
ans =
15 34
18 32
43 74

(3) Matrix Inversion

The function inv() is used to calculate the inverse of the matrix A. Note that, the matrix A must be a square
matrix (same number of rows and columns) to have an inverse matrix.
Example 10
>> A = [2 1 4; 4 1 8; 2 -1 3];
>> inv(A)
ans =
5.5000 -3.5000 2.0000
2.0000 -1.0000 0
-3.0000 2.0000 -1.0000
Note: The Command A^(-1)will have the same effect.

(4) Matrix Division: Left Division and Right Division

There are two kinds of matrix division in MATLAB for two matrices A and B – Left division and Right
division. Mathematically, for two matrices A and B, left division is the operation: and right division is
the operation: .

10
 For illustration, let,

[ ] [ ]

In MATLAB –
>> A = [1 2 3; 4 5 6; 7 8 9];
>> B = [10 11 12; 13 14 15; 16 17 18];

Mathematical
Command Name Example
Operation
>> A\B
ans =
A\B Left Division -20.0000 -17.0000 -18.0000
33.0000 26.0000 27.0000
-12.0000 -8.0000 -8.0000
>> A/B
ans =
A/B Right Division 1.0000 3.0000 -3.0000
1.0000 2.0000 -2.0000
0.5000 2.0000 -1.5000

(5) Matrix Element-by-Element Operation (Very Very Important!)

Let, two matrices A and B are defined as –

[ ] [ ]

Then, the element-by-element operations are:

Command Meaning Outcome

Multiplies every element of A with the
A.*B [ ]
corresponding element of B

Divides every element of A by the

A./B [ ]
corresponding element of B

A.^n Every element of A is raised by a power n [ ]

Note: Observe carefully the difference between using dot and not using dot.

(6) Built-in-functions for Matrix Analysis

Function Description Example

>> A = [5 9 2 4];
>> mean(A)
mean(A) Returns mean value of the elements of vector A
ans =
5

11
 >> A = [5 9 2 4];
Returns the median value of the >> median(A)
median(A)
elements of vector A ans =
4.5000
>> A = [5 9 2 4];
Returns maximum value of the >> a = max(A)
a = max(A)
elements of vector A a =
9
>> [a, m] = max(A)
a =
Returns maximum value of the elements of 9
[a, m] = max(A)
vector A and also returns its position m =
2
>> A = [5 9 2 4];
Returns minimum value of the >> b = min(A)
b = min(A)
elements of vector A a =
2
>> [b, n] = min(A)
b =
Returns minimum value of the elements of 2
[b, n] = min(A)
vector A and also return its position n =
3
>> A = [5 9 2 4];
>> sum(A)
sum(A) Returns the sum of the vector A
ans =
20
>> A = [5 9 2 4];
Arranges the elements of vector A >> sort(A)
sort(A)
in their ascending order ans =
2 4 5 9
>> A = [5 9 2 4];
Returns the standard deviation of >> std(A)
std(A)
the elements of vector A ans =
2.9439
>>A = [2 4 6;…
6 9 1;…
Returns the determinant value of matrix A. 0 2 7];
det(A)
Here, A must be a square matrix >> det(A)
ans =
26
>> A = [1 2 3];
Returns the scalar or dot product of >> B = [3 4 5];
dot(A, B)
two vectors A and B >> dot(A, B)
ans =
26
Returns the vector or cross product of >> cross(A, B)
cross(A, B) ans =
two vectors A and B
-2 4 -2

12
REPORT
[Submit as in the next session]
Q1. Write a MATLAB program to convert the given temperature in Degree Celsius to Degree Fahrenheit.
[ ]

Q2. Write a MATLAB program to evaluate the efficiency, :

Where,

Q3. For an RLC circuit, the damped natural frequency is given by:

Where, , and . Write a program to calculate the for the values of varying from
in a step of .

Q4. Use matrix operations to solve the following system of linear equations:

[Hint: where, [ ] and [ ]. Use Matrix division then.]

Reference Books:
(1) “Duane Hanselman, Bruce Littlefield “, Mastering MATLAB 7
(2) “Amos Gilat”, MATLAB An introduction with Applications
Web resource: [Link]

13
 Daffodil International University
Department of Electrical and Electronic Engineering
EEE 322: Digital Signal Processing Laboratory

EXPERIMENT NO: 01-B

NAME OF THE EXPERIMENT: SCRIPT FILES IN MATLAB

Why Script Files

In the previous experiment, the commands were executed in the command window. Although every MATLAB
command can be executed in this way, there are some drawbacks of using command window:
(i) Using the command window to execute a series of commands (especially if they are related to each other) is not
convenient and may be difficult or even impossible.
(ii) The commands in the command window cannot be saved and executed again.
(iii) The command window is not interactive. This means that, every time the ENTER key is pressed, only the last
command is executed and everything executed before is unchanged. If a change or correction is needed in a
previous command,the commands have to be entered again.
A different way of executing commands with MATLAB is first to create a file with a list of commands, save it and
then run the file. The files that are used for this purpose is called script files.

Creating and Saving a script file

In MATLAB, script files are created and edited in the Editor/Debugger window. To do so, follow: File New M-
[Link] will see the following window:

|
Run icon

The commands in the script file are

typed line by [Link] lines are
numbered automatically.A new line
Line number starts when Enter key is pressed.

Figure 2.1: The editor/debugger window

14
Before a script file can be executed it has to be saved. This is done by choosing save as from the file menu, selecting
a location and entering a name for the file.

Running a script file

A script file can be executed in several ways:
Directly from the editor window –
(i) Click the Run icon (see Fig. 2.1)
(ii) Follow: debug  run <file_name.m>
(iii) Press F5 button in keyboard
Or, indirectly from the command window –
Type filename in the command window and press Enter

Example 1
Voltage Divider: In an electrical circuit if several resistors are connected in series, the voltage across each of them
is given by the voltage divider rule:

R1 R2 R3
+
vs _
R4
Where, and are the voltage across nth resistor, .
R7 R6 R5
is the equivalent resistance of the connection ( ∑ ).
The power dissipated in the nth resistor is given by: Figure 2.2: A Series Circuit

4 5

Write a program in script file, which calculates the voltage across each resistor and power dissipated in each resistor,
in a series circuit in Fig. 2.2. When the script file is executed, it prompts the user to first enter the source voltage and
then to enter the resistance of the resistors in a vector. The program displays the result in a table. Following the
table, the program displays the current in the circuit, and the total power. Run the program and enter the following
data for and the ‟s.

Solution: Create the m-file as in the editor window given below:

15

Figure 2.3: Example of a script file

After executing the script file in the command window ,the output is as follows:

>> Vdivider
Please enter the source voltage:24
Please enter the values of the resistors in a row vector
[20 14 12 18 8 15 10]
Resistance Voltage Power
(Ohms) (Volts) (Watts)
20.0000 4.9485 1.2244
14.0000 3.4639 0.8571
12.0000 2.9691 0.7346
18.0000 4.4536 1.1019
8.0000 1.9794 0.4897
15.0000 3.7113 0.9183
10.0000 2.4742 0.6122
The current in the circuit is 0.247423 Amps.
The total power dissipated in the circuit is 5.938144 Watts.

Controlling Execution Flow

It is possible to do a great deal in MATLAB by simply executing statements involving matrix expressions, one after
the other. However, there are cases in which one simply must substitute some non-sequential order. To facilitate
this, MATLAB provides three relatively standard methods for controlling program flow: the if-else statement,
the for loop and the while loop.

Examples of control flow statements

if-else
num = -3;
if num > 0
x = ’positive’;
elseif num < 0
x = ’negative’;
elseif num == 0
x = ’zero’;
else
x = ’error’;
end
Question: What is the value of ‟x‟ after execution of the above code?

while
x = -10;
while x < 0
x = x + 1;
end
Question: What is the value of x after execution of the above code?

for
x = 0;
for i = 1:10
x = x + 1;
end

16
Question: What is the value of x after executing the above code?

switch
x = 10;
switch x
case 5
disp(‘x is 5’)
case 6
disp(‘x is 6’)
case 8
disp(‘x is 8’)
case 10
disp(‘x is 10’)
otherwise
disp(‘invalid choice’)
end
Question: What is the output shown after executing the above code?

break
The break statement lets you exit early from a for or while loop:
x = -10;
while x < 0
x = x + 2;
if x == -2
break;
end
end
Question: What is the value of x after executing the above code?

Relational and Logical Operators

MATLAB supports the following relational and logical operators:
Relational Operators Logical Operators
Symbol Meaning Symbol Meaning
<= Less than or equal to & AND
< Less than | OR
>= Greater than or equal to ~ NOT (unlike C. In C, it was!‟)
> Greater than
== Equal to
Not equal to
~=
(unlike C. In C, it was „!=‟)

17
Example 2
Write a script file to sum the series , taken from user.

Solution:

Figure 2.4: Use of for loop

Executing from command window

Figure 2.5: Execution of series.m

18
Writing User-defined Functions
A function file starts with a line declaring the function with its name, input arguments and outputs. Then, follows the
statements required to produce the outputs from the inputs.
Example: Suppose, we want to calculate average of some numbers. For this purpose, let us create a function named
average.m. The steps are:
(i) Open a new m-file from editor.
(ii) Write the following code as shown in the figure:

Must be saved as the same name of

the function

Figure 2.6: Example of a function script file.

Structure of a function definition line

function [output_arguments] = function_name(input_arguments)

Every function definition line List of output arguments

must start with the keyword typed inside parentheses
function (all lower case)

List of input arguments

typed inside parentheses
The name of the function
This name is very important as the function will be
called by this name. Also the script file has to be
saved by the same name i.e., function_name.m

Body of a function
The function body contains the computer program (code) that actually performs the computation. The code can use
all MATLAB programming features. This includes calculations, assignments, any built-in or another user-defined
functions, flow control etc. You can now understand the function as shown in Fig. 2.6.

19
How to use a user-defined function
We can call a function from the command window as show in the following Fig. 2.7.

Figure 2.7: Calling our function average.m from command window

Example 3
Write a function that takes radius as input and find the area and perimeter of a circle.
Solution: The function is defined as follows:

function [area, perimeter] = circle(radius)

area = pi*radius^2;
perimeter = 2*pi*radius;

Create it in m-file script and save as circle.m

Calling circle.m function from command window:

>> r = 10;
>> [A, C] = circle(r)
A =
314.1593
C =
62.8319

2-D Plotting in MATLAB

Plots are very useful tools for presenting information. MATLAB has many commands that can be used for creating
different types of plots. We will discuss some of them now.

20
plot command
The plot command is used to create two-dimensional plots. The simplest form of the command is:
plot(x,y)
Example: plot command to draw simple curve.
>> x = [1, 2, 3, 5, 7, 8, 10];
>> y = [2, -6.5, 7, 7, 5.5, 4, 6];
>> plot(x, y)

The commands on the command window give the

following figure:

The plot command has additional optional arguments that

can be used to specify the color and style of the line and the
color and type of markers,if any are [Link] these
options the command has the form:

plot(x, y, ’LineSpecifiers’, ’PropertyName’, ’PropertyValue’)

Line Specifiers
Line specifiers are optional and can be used to define the style and color of the line and type of the markers (if the
markers are used).
Line-style Specifiers Line-color Specifiers
Used to define the style of the line i.e. – solid line, Used to specify the color of the line:
dashed line and so on.
Line-color Specifier
Line-style Specifier
Red r
Solid (default) - g
Green
Dashed -- b
Blue
Dotted . c
Cyan
Dash-dot -. m
Magenta
Pointed o y
Yellow
Black k
White W
Marker-type Specifiers
Marker type specifiers are used to define the type of the markers (when markers are used).
Marker-type Specifier
Plus Sign +
Circle o
Asterisk *
Point .
Square s
Diamond d
Five-point Star p
Six-point Star H

21
 Some Examples

Command Meaning
plot(x, y) A blue solid line connecting the points with no markers (default)
plot(x, y, ’r’) A red solid line connecting the points.
plot(x, y, ’--y’) A yellow dashed line connecting the lines
plot(x, y, ’*’) The points are marked with * (no line between points)
A green dotted line connecting the points that are marked with
plot(x, y, ’g:d’)
diamond markers.

Example 4
Suppose the following given data is used to create vectors that are then used in the plot command.
Year 2002 2003 2004 2005 2006 2007 2008
Sales(millions) 8 12 20 22 18 24 27
To plot this data, the following codes are written in the
command window:

>> yr = 2002:2008;
>> sle = [8 12 20 22 18 24 27];
>> plot(yr, sle, '--r*',...
'linewidth', 2,...
'markersize', 12)

And the output figure will be as the following:

Plot of Given Functions

In many situations, there is a need to plot a given function.
This can be done in MATLAB by using plot or fplot function.
(i) Using plot function
Suppose we want to plot the function –
( ) ( )
The program that plots this function is shown in the
following commands:
>> x = -2:.001:4;
>> y = 3.5.^(-0.5*x).*cos(6*x);
>> plot(x, y)

And the output figure will be as the following:

(ii) Using fplot function

The commands will be:
>> fplot('3.5.^(-0.5*x).*cos(6*x)',[-2 4])

The output figure is the same as above.

22
Plotting Multiple Graphs in the same plot
There are three ways to do so:
(i) Using plot command with pair of vectors inside it
(ii) Using hold on and hold off command
(iii) Using line command

Example 5
Plot the function and its first and second derivatives for , all in the same plot.

Solution:
x = -2:0.01:4;
y = 3*x.^3-26*x+6;
yd = 9*x.^2-26;
ydd = 18*x;

% using ‘plot’ command

plot(x, y, '-b',...
x, yd, '--r',...
x, ydd, ':k')

% using ‘hold’ command

figure(2)
plot(x, y, '-b')
hold on
plot(x, yd, '--r')
plot(x, ydd, ':k')
hold off

% using ‘line’ command Figure 2.8: Plotting multiple functions in a single graph
figure(3)
plot(x, y, 'linestyle', '-', 'color', 'b')
line(x, yd, 'linestyle', '--', 'color', 'r')
Any one of the three commands will give the desiredplot in a single figure as shown in Fig. 2.8.

Formatting a plot using commands

The formatting commands are placed after the plot or fplot commands .The various formatting commands are
shown in the following table:
Labeling the Axes
Command Purpose
xlabel(‘string’)
To label the x-axis and y-axis by text.
ylabel(‘string’)
title(‘string’) To give a title to your figure
text(x, y, ’string’)
To place a text anywhere on the figure.
gtext(‘string’)
legend(‘string1’, ’string2’) To place a legend on the plot

23
 Setting Limits on Axes
Command Purpose
Sets the limits of the x-axis
axis([xmin, xmax])
(xmin and xmax are numbers)
axis([xmin, xmax, ymin, ymax]) Sets the limits of both the axis
axis equal Sets the same scale for both the axis
axis square Sets the axis region to be square

Setting Gridline ON/OFF

Command Purpose
grid on Adds grid line to the plot
grid off Removes grid line from the plot

Making plot of discrete-time plots

Use of stem
The function stem is used to plot discrete-time
functions (in contrast to the plot function, which is
used to plot continuous-time functions).
Example 6
Observe the following code:
>> x = [Link];
>> y = sin(x);
>> stem(x, y)

This will give the plot in Fig. 2.9.

Figure 2.9: Plotting discrete-time functions

Plotting Multiple Plots on the Same Figure

Use of subplot
Multiple plots on the same page can be created with
(3,2,1) (3,2,2)
the subplot command which has the form:
subplot(m, n, k)

The command divides the figure window (page when (3,2,3) (3,2,4)
printed) into m × n rectangular subplots where plots
will be created. k varies from 1 to mn row-wise.
(3,2,5) (3,2,6)
For example, the commands subplot(3,2,1) …
subplot(3,2,6)will create the arrangement as
shown in Fig. 2.10. Figure 2.10: Use of subplot to divide the
figure window into multiple subfigures

24
Example 7
The figure shows a simple high-pass filter consisting of a resistor and a capacitor. The ratio of the output voltage
to the input is given by:
+

+
C
Vi R Vo
_
Assume and . Calculate and plot the amplitude and _

phase response of this filter as a function of frequency.

Solution:The script file is as follows:

The output figure is as follows:

25
REPORT
[Submit in the next session]
Q1. Write MATLAB script to plot and its derivative on the same plot in two different colors and in three
different ways.

Q2. A resistor of value and an inductor of value are connected in series to a voltage source, as
shown in the figure (a). V

R 12V

+
vs _
L

0 t
0.5s
(a) (b)

When the voltage source applies a rectangular voltage pulse with an amplitude of and a duration of
, as shown in figure (b). The current ( ) is given by:
. /
( ) {
. /
Make a plot of the current as a function of time for .

Q3. The figure shows a simple low-pass filter consisting of a resistor and a capacitor.

R +

+ V0
vi C
_
_

Assume and , where is the last two digits of your student ID. Calculate the amplitude and
phase response of this filter as a function of frequency and plot them in a single figure using subplot.
[Hint:Gain of a low-pass filter is given by: ].

Reference Books:
(1) “Duane Hanselman, Bruce Littlefield “, Mastering MATLAB 7
(2) “Amos Gilat”, MATLAB An introduction with Applications
Web resource: [Link]

26
 Daffodil International University
Department of Electrical and Electronic Engineering
EEE 322: Digital Signal Processing Laboratory

EXPERIMENT NO: 02

NAME OF THE EXPERIMENT: IMPLEMENTATION OF DISCRETE-TIME SIGNALS

Generation of DT signals
In MATLAB, we can represent a finite-duration sequence by a row vector of appropriate values. However, such a
vector does not have any information about sample position . Therefore, a correct representation of ( ) would
require two vectors one each for and .
For example, a sequence ( ) 2 3 can be represented in MATLAB as:
>> n = [-3, -2, -1, 0, 1, 2, 3, 4];
>> x = [2, 1, -1, 0, 1, 4, 3, 7];
>> stem(n,x)
The corresponding DT signal is shown in the following figure:

Generation of Common Sequences

A number of basic sequences are commonly used in Discrete Signal Processing applications. They are discussed in
the following table:
Function Mathematical
SL MATLAB Function
Name Definition
function [x, n] = delta(n0, n1, n2)
Unit Sample
1. ( ) 2 n = n1:n2;
Function x = [(n - n0) == 0];
function [x, n] = step(n0, n1, n2)
Unit Step
2. ( ) 2 n = n1:n2;
Function x = [(n - n0) >= 0];

27
 Real valued function [x, n] = realx(a, n1, n2)
( )
3. Exponential n = n1:n2;
Sequence x = a.^n;
function[y, n] = sigadd(x1, n1, x2, n2)
n_start = min(min(n1), min(n2));
n_end = max(max(n1), max(n2));
n = n_start:n_end;
y1 = zeros(1, length(n));
Signal * ( )+ * ( )+
4. y2 = y1;
Addition * ( ) ( )+
n11 = find((n>=min(n1))&(n<=max(n1))==1);
n22 = find((n>=min(n2))&(n<=max(n2))==1);
y1(n11) = x1;
y2(n22) = x2;
y = y1 + y2;

function [y, n] = sigmult(x1, n1, x2, n2)

n_start = min(min(n1), min(n2));
n_end = max(max(n1), max(n2));
n = n_start:n_end;
y1 = zeros(1, length(n));
Signal * ( )+ * ( )+
5. y2 = y1;
Multiplication * ( ) ( )+
n11 = find((n>=min(n1))&(n<=max(n1))==1);
n22 = find((n>=min(n2))&(n<=max(n2))==1);
y1(n11) = x1;
y2(n22) = x2;
y = y1.*y2;
function [y, n] = sigshift(x, m, n0)
6. Signal Shifting ( ) * ( )+ n = m + n0;
y = x;

function [y, n] = sigfold(x, n)

7. Signal Folding ( ) * ( )+ y = fliplr(x);
n = -fliplr(n);

Note: fliplr is a library function. Write help fliplr in command window to know in details.

Illustrative Examples
Example 1
Generate and plot the following sequences over the indicated intervals:
(a) ( ) ( ) ( )
( )
(b) ( ) , ( ) ( )- , ( ) ( )-
(c) ( ) ( ) ( )
where, ( ) is the Gaussian random sequencewith zero mean and unit variance .

Solution: The necessary code is listed below:

% for part a
n = -5:5;
x = 2*delta(-2,-5,5) - delta(4,-5,5);
subplot(3,1,1), stem(n,x)
title('sequence in part a.')
xlabel('n'), ylabel('x(n)')

28
% for part b
n = 0:20;
x1 = n.*(step(0,0,20) - step(10,0,20));
x2 = 10*exp(-0.3*(n - 10)).*(step(10,0,20) - step(20,0,20));
x = x1 + x2;
subplot(3,1,2), stem(n,x)
title('sequence in part b.')
xlabel('n'),ylabel('x(n)')

% for part c
n = 0:50;
w = randn(size(n));
x = cos(0.04*pi*n) + 0.2*w;
subplot(3,1,3), stem(n,x)
title('sequence in part c.')
xlabel('n'),ylabel('x(n)')

The output is shown below:

sequence in part a.
2
x(n)

-2
-5 -4 -3 -2 -1 0 1 2 3 4 5
n
sequence in part b.
10
x(n)

0
0 2 4 6 8 10 12 14 16 18 20
n
sequence in part c.
2
x(n)

-2
0 5 10 15 20 25 30 35 40 45 50
n

29
Example 2
Generate and plot the samples (use the stem function) of the following sequence using MATLAB:

(a) ( ) ∑( ), ( ) ( )-

(b) ( ) , ( ) ( )- ( ) ( ) , ( ) ( )-
(c) ( ) ( ) . /
(d) ( ) ( ) ( )
Where, ( ) is a random sequence uniformly distributed between , -
(e) ( ) * + plot periods.

Solution: The necessary code is listed below:

% (a)
n = 0:25;
x1 = zeros(1,26);
for m=0:10,
x1 = x1 + (m+1).*(delta(2*m,0,25)- delta(2*m+1,0,25));
end
subplot(5,1,1), stem(n, x1)
title('Sequence in Problem a')

% (b)
n = -25:25;
x2 = n.^2.*(step(-5,-25,25)- step(6,-25,25)) + 10*delta(0,-25,25) +...
20*0.5.^n.*(step(4,-25,25)- step(10,-25,25));
subplot(5,1,2), stem(n, x2)
title('Sequence in Problem b')

% (c)
n = 0:20;
x3 = 0.9.^n.*cos(0.2*pi.*n + pi/3);
subplot(5,1,3), stem(n,x3)
title('Sequence in Problem c')

% (d)
n = 0:100;
w = randn(size(n));
x4 = 10.*cos(0.008*pi.*n.^2) + w;
subplot(5,1,4), stem(n, x4)
title('Sequence in Problem d')

% (e)
x5 = [2,1,2,3];
x5 = x5'*ones(1,5); % Plots 5 periods
x5 = (x5(:))';
n = 0:size(x5);
subplot(5,1,5), stem(n, x5)
title('Sequence in Problem e')

30
The corresponding output is shown in the following figure:

Sequence in Problem a
10
0
-10
0 5 10 15 20 25
Sequence in Problem b
40
20
0
-25 -20 -15 -10 -5 0 5 10 15 20 25
Sequence in Problem c
1
0
-1
0 2 4 6 8 10 12 14 16 18 20
Sequence in Problem d
20
0
-20
0 10 20 30 40 50 60 70 80 90 100
Sequence in Problem e
4
2
0
0 2 4 6 8 10 12 14 16 18 20

Example 3
Let ( ) 2 3
Generate and plot the following sequence using MATLAB:
(a) ( ) ( ) ( )
(b) ( ) ( ) ( ) ( )

Solution: The necessary code is listed below:

n = -2:10;
x = [1:7, 6:(-1):1];
% (a) x1(n) = 2*x(n-5) - 3*x(n+4)
[x11,n11] = sigshift(x,n,5);
[x12,n12] = sigshift(x,n,-4);
[x1,n1] = sigadd(2*x11,n11,-3*x12,n12);
subplot(2,1,1), stem(n1, x1)
title('Sequence in a')
xlabel('n'), ylabel('x1(n)')
axis([min(n1)-1, max(n1)+1, min(x1)-1, max(x1)+1])

% (b) x2(n) = x(3-n) + x(n)*x(n-2)

[x21,n21] = sigfold(x,n);
[x21,n21] = sigshift(x21,n21,3);

31
 [x22,n22] = sigshift(x,n,2);
[x22,n22] = sigmult(x,n,x22,n22);
[x2,n2] = sigadd(x21,n21,x22,n22);
subplot(2,1,2), stem(n2, x2)
title('Sequence in Example b')
xlabel('n'), ylabel('x2(n)')
axis([min(n2)-1, max(n2)+1, 0, 40])

The corresponding output is:

Sequence in a

10

0
x1(n)

-10

-20
-6 -4 -2 0 2 4 6 8 10 12 14 16
n
Sequence in b
40

30
x2(n)

20

10

0
-8 -6 -4 -2 0 2 4 6 8 10 12
n

Example 4
( )
Generate the complex valued signal, ( ) and plot its magnitude,phase,the real part
and the imaginary part in four separate subplots.

Solution: The necessary code is listed below:

n = -[Link];
alpha = -0.1 + 0.3j;
x = exp(alpha*n);

subplot(2,2,1), stem(n, real(x))

title('real part'), xlabel('n')
subplot(2,2,2), stem(n, imag(x))
title('imaginary part'), xlabel('n')
subplot(2,2,3), stem(n, abs(x))
title('magnitude part'), xlabel('n')
subplot(2,2,4), stem(n,(180/pi)*angle(x))
title('phase part'), xlabel('n')

32
The corresponding output is:

real part imaginary part

2 1

0 0

-2 -1

-4 -2
-10 -5 0 5 10 -10 -5 0 5 10
n n
magnitude part phase part
3 200

100
2
0
1
-100

0 -200
-10 -5 0 5 10 -10 -5 0 5 10
n n

REPORT
[Submit on the next lab session]
Q1. A unit ramp DT sequencer(n) is defined as:
( ) 2
Develop a function of the form: function [x,n] = r(n)to implement the above [Link] the above
sequence to plot a ramp sequence for the interval to .

Q2. Construct a function which decomposes a real signal into even and odd parts having the following heading:
function [xe, xo, m] = evenodd(x, n)
where, xe and xo are the even and odd part of the sequence x and m = m1:m2, m1 & m2 are the start point of
and end point, respectively, of the even/odd part. Use this function to decompose ( ) ( ) ( ).
Show original sequence, even part and odd part in subplots.

Q3. Generate and plot the samples (use the stem function) of the following sequence using MATLAB:
( ) ( ) ( ) ( )where, ( ) 2 3

Reference Books:
(1) “S.K. Mitra “, Digital Signal Processing:A computer based Approach, 3rd Edition.
(2) “J.G. Proakis and D.G. Manolakis”, Digital Signal Processing: Principles, Algorithms and Applications,
4th Edition
Web resource: [Link]

33
 Daffodil International University
Department of Electrical and Electronic Engineering
EEE 322: Digital Signal Processing Laboratory

EXPERIMENT NO: 03

NAME OF THE EXPERIMENT: SAMPLING & RECONSTRUCTION OF ANALOG

SIGNALS

Introduction
Signals available in our environment are mostly analog in nature. As a result, we cannot use them for digital signal
processing as DSP deals with discrete-time signals. As DSP is unable to handle continuous-time signals, we need to
use sampling to convert the continuous-time signals into discrete-time signals for processing.

Theory
Sampling is the process of converting a continuous-time signal into the corresponding discrete-time signal. For the
sake of conversion, samples of the continuous-time signal are taken periodically. The interval between two sampling
instant is called the sampling period.
If ( ) is a continuous time signal, the corresponding discrete signal will be ( ) where, and is
the sampling period. We can also write the discrete signal as ( ).
( ) ( ) ( )

Figure 4.1: Illustration of Sampling of an Analog Signal

Sampling Rate
Sampling rate (or, sampling frequency) is the number of samples taken in a second by the sampler. If the sampling
period is , then the sampling rate will be:

( )
If the sampling rate is increased, the number of samples taken per second will be increased and vice versa.
Hence, in the case of a sinusoidal signal, ( ) ( ):
( ) ( ) , ( ) -

[ ( ) ] ( )
( )
Where, is the digital frequency corresponding to the analog frequency, .

34
Sampling Theorem
The exact reconstruction of a continuous-time baseband signal from its samples is possible if the signal is band-
limited and the sampling rate is greater than twice the signal bandwidth.
Let, the signal ( ) has a maximum frequency of , i.e., the signal bandwidth, . Then, the sampling
theorem dictates that, for successful reconstruction of ( ) from its sample ( ) –
( )
The sampling rate is regarded as the Nyquist rate, .
This theorem also leads to an effective reconstruction formula:

̂ ( ) ∑ ( ) ( )
( )
( )
( ) ( )
Where, ( ) is the ideal interpolation function as illustrated in Fig. 4.2.

Figure 4.2: The interpolation function— sinc function

Since, . / ( ) ( ), then the reconstruction formula in ( ) can be rewritten as –

̂ ( ) ∑ ( ) , ( )- ∑ ( ) , ( )- ( )

Here, the formulas in ( ) and ( ) are used in the examples to perform sampling and reconstruction, respectively.

MATLAB Examples
Sampling & Reconstruction of a Sinusoid
Example 1
Consider an analog signal, ( ) ( ) . It is sampled at
intervals to obtain ( ). Use .
(a) For each , plot ( ).
(b) Reconstruct the analog signal ( ) from the samples ( ) by means of of the sinc interpolation. Estimate the
frequency in ( ) from your plot.
(c) Comment on your results.

35
Solution:
(a) MATLAB Code for Sampling
% Analog signal...
t = 0:0.001:1;
xa = cos(20*pi*t); % F = 10 Hz
subplot(4,1,1), plot(t, xa, 'k')
ylabel('x_a(t)'), title('Analog Signal, x_a(t) = cos(20\pit)')

% Sampling...
T = 0.01; N1 = round(1/T);
n1 = 0:N1;
x1 = cos(20*pi*n1*T); % x(n) = xa(nT)
subplot(4,1,2), stem(n1*T, x1, 'r')
axis([0, 1, -1.1, 1.1]); ylabel('x_1(n)')
title('Sampling of x_a(t) using T = 0.01 sec')

T = 0.05; N2 = round(1/T);
n2 = 0:N2;
x2 = cos(20*pi*n2*T);
subplot(4,1,3), stem(n2*T, x2, 'm')
axis([0, 1, -1.1, 1.1]); ylabel('x_2(n)')
title('Sampling of x_a(t) using T = 0.05 sec')

T = 0.1; N3 = round(1/T);
n3 = 0:N3;
x3 = cos(20*pi*n3*T);
subplot(4,1,4), stem(n3*T, x3, 'b')
axis([0, 1, -1.1, 1.1]); ylabel('x_3(n)')
title('Sampling of x_a(t) using T = 0.1 sec')

(b) MATLAB Code for Reconstruction

% Reconstruction...
figure(2)
subplot(4,1,1), plot(t, xa, 'k')
ylabel('x_a(t)'), title('Original Analog Signal, x_a(t) = cos(20\pit)')

T = 0.01; Fs = 1/T;
ya1 = x1*sinc(Fs*(ones(length(n1), 1)*t - (n1*T)'*ones(1, length(t))));
subplot(4,1,2), plot(t, ya1, 'r')
axis([0, 1, -1.1, 1.1]); ylabel('y_a(t)')
title('Reconstruction of x_a(t) when T = 0.01 sec')

T = 0.05; Fs = 1/T;
ya2 = x2*sinc(Fs*(ones(length(n2), 1)*t - (n2*T)'*ones(1, length(t))));
subplot(4,1,3), plot(t, ya2, 'm')
axis([0, 1, -1.1, 1.1]); ylabel('y_a(t)')
title('Reconstruction of x_a(t) when T = 0.05 sec')

T = 0.1; Fs = 1/T;
ya3 = x3*sinc(Fs*(ones(length(n3), 1)*t - (n3*T)'*ones(1, length(t))));
subplot(4,1,4), plot(t, ya3, 'b')
axis([0, 1, -1.1, 1.1]); ylabel('y_a(t)')
title('Reconstruction of x_a(t) when T = 0.1 sec')

36
 The output figure for part (a) is shown below:

The output figure for part (b) is shown below:

(c) Comments: From the plots, it is clear that reconstructions from samples at depict
the original frequency, but the reconstruction for shows the original frequency aliased to zero.

37
Observing the Effect of Aliasing
Example 2
Consider the two following signals for :
(i) ( )
(ii) ( )
(a) Plot ( ) and ( ) on the same plot. Use .
(b) Sample both signals at to obtain ( )and ( ). Plot these two signals on the same plot also.
(c) What does the result indicate? Explain.

Solution:
MATLAB Code
t = -0.5:0.001:0.5;
xa1 = cos(6*pi*t);
xa2 = cos(14*pi*t);

T = 0.1; N = round(1/T); % Fs = 10 Hz
n = -N/2:N/2;
x1 = cos(6*pi*n*T); % x(n) = xa(nT)
x2 = cos(14*pi*n*T);

subplot(2,1,1), plot(t, xa1, 'k'), hold on

stem(n*T, x1, 'b'), axis([-0.5, 0.5, -1.1, 1.1])
ylabel('Amplitude'), title('Sampling of x_{a_1}(t) using T = 0.1 sec')
legend('x_{a_1}(t)', 'x_1(n)')
subplot(2,1,2), plot(t, xa2, 'k'), hold on
stem(n*T, x2, 'b'), axis([-0.5, 0.5, -1.1, 1.1])
ylabel('Amplitude '), title('Sampling of x_{a_2}(t) using T = 0.1 sec')
legend('x_{a_2}(t)', 'x_2(n)')

The resultant plot for part (a) and (b) is:

Note that, even though the two signals, ( ) and ( ) are distinct in time-domain, the sample points obtained
from them are undistinguishable (For or, ).

38
REPORT
[Submit on the next session]
Q1. An analog signal ( ) ( )is sampled using the following sampling periods:
(a)
(b)
(c)
In each case, plot the spectrum in the resulting discrete-time signal.

| |
Q2. Let, ( ) . Sample the function at and comment on
your result.

Reference Books:
(1) “S.K. Mitra”, Digital Signal Processing: A computer based Approach, 3rd Edition.
(2) “J.G. Proakis and D.G. Manolakis”, Digital Signal Processing: Principles, Algorithms and Applications,
4th Edition
Web resource: [Link]

39
 Daffodil International University
Department of Electrical and Electronic Engineering
EEE 322: Digital Signal Processing Laboratory

EXPERIMENT NO: 04

NAME OF THE EXPERIMENT: QUANTIZATION OF DISCRETE-TIME SIGNALS

Introduction
The basic task of the A/D converter is to convert a continuous range of input amplitudes into a discrete set of digital
code-words. Along with sampling, this conversion also involves the processes of quantization and coding. As DSP is
unable to handle analog signals (continuous-time and continuous-valued signal), the sampling process is applied to
discretize the time-axis. However, the amplitudes of the resultant signals still have continuous ranges of values.
Quantization is used to discretize the amplitude axis and process them using DSP. It converts a continuous-
amplitude sample to a discrete-amplitude sample.

Theory
In digital signal processing, quantization is the process of approximating a continuous range of values by a
relatively small set of discrete values. Quantization is a nonlinear and noninvertible process that maps a given
amplitude ( ) ( ) at time into an amplitude , taken from a finite set of values. All sample values
falling between two quantization levels are considered to be located at the center of the quantization interval. The
procedure is illustrated in Fig. 5.1(a), where the signal amplitude range is divided into intervals.
* ( ) + ( )
by decision levels .

Figure 5.1: (a) Quantization Process and (b) A Mid-tread Quantizer

The possible outputs of the quantizer (i.e., the quantization levels) are denoted as: ̂ ̂ ̂ . The operation of the
quantizer is defined by the relation:
( ) * ( )+ ̂ ( ) ( )

40
Uniform Quantizer
The most common and most used form of quantization is the uniform or linear quantization, defined by
( )
Where, is the quantizer step-size.
For a discrete-time signal, ( ) with a maximum and minimum amplitudes of and , the quantization step-
size (for levels) is defined as:

( )
Where, is the dynamic range of the uniform quantizer.
The input-output characteristics of a uniform quantizer are shown in the following figure:

Figure 5.2: The I/O Characteristics of an Uniform Quantizer

In quantization, all sample values falling between two quantization levels are considered to be located at the center
of the quantization interval. In this manner, the quantization process introduces a certain amount of error or
distortion into the signal samples, known as the quantization error or,quantization noise.
( ) ( ) ( ) ( )
The maximum magnitude of quantization noise is: , i.e., ( ) . This noise is minimized by establishing
a large number of small quantization intervals. As all the quantization levels are assigned to distinct binary code-
words, the increased number of quantization levels also increases the number of bits for a single code-word. For
example, for a number of levels, , the number of bits is .

MATLAB Examples
Quantization of vectors
Example 1
(a) Write a function uquant to quantize a vector, in the form:
function y = uquant(x, L)
where is the number of levels..
(b) Quantize the vector: ( ) , - for the following number of
levels, .

Solution:

41
(a) MATLAB Code for uquant Function
function y = uquant(x, L)
xmax = max(max(x));
xmin = min(min(x));
stepsize = (xmax - xmin)/(L - 1);
y = round((x - xmin)/stepsize)*stepsize + xmin;

(b) MATLAB Code for Quantization

x = [1 2 -3 -4 5 10 12 13 16 11 0 -5 7 9 10 16]
z = length(x);
n = 1:z;

y1 = uquant(x, 2)
subplot(2,2,1), stem(n, y1)
axis([1, z, min(x)-1, max(x)+1]), ylabel('x_q(n)')
title('No. of levels, L = 2')

y2 = uquant(x, 4)
subplot(2,2,2), stem(n, y2)
axis([1, z, min(x)-1, max(x)+1]), ylabel('x_q(n)')
title('No. of levels, L = 4')

y3 = uquant(x, 8)
subplot(2,2,3), stem(n, y3)
axis([1, z, min(x)-1, max(x)+1]), ylabel('x_q(n)')
title('No. of levels, L = 8')

y4 = uquant(x, 16)
subplot(2,2,4), stem(n, y4)
axis([1, z, min(x)-1, max(x)+1]), ylabel('x_q(n)')
title('No. of levels, L = 16')

The data output is:

x =
1 2 -3 -4 5 10 12 13 16 11 0 -5 7 9 10 16

y1 =
-5 -5 -5 -5 -5 16 16 16 16 16 -5 -5 16 16 16 16

y2 =
2 2 -5 -5 2 9 9 16 16 9 2 -5 9 9 9 16

y3 =
1 1 -2 -5 4 10 13 13 16 10 1 -5 7 10 10 16

y4 =
Columns 1 through 8
0.6000 2.0000 -3.6000 -3.6000 4.8000 10.4000 11.8000 13.2000

Columns 9 through 16
16.0000 10.4000 0.6000 -5.0000 7.6000 9.0000 10.4000 16.0000

42
The figure output is:

Quantization of sinusoidal signals

Quantization of sinusoidal signals is similar to the quantization of a vector. Here, the vector again contains the
amplitude of the signal at various time instants.

Example 2
Quantize the signal, ( ) ( ) for the following number of levels, and
compare with the original analog signal in the same plot. Take points.

Solution:
N = 21;
t = linspace(-1, 1, N);
Y = 21*cos(2*pi*t);

Y1 = uquant(Y, 4);
subplot(2,2,1), plot(t, Y, 'k'), hold on
stem(t, Y1), ylabel('Amplitude'), title('No. of levels, L = 4')
legend('Y_a(t)', 'Y_q[n]')

Y2 = uquant(Y, 8);
subplot(2,2,2), plot(t, Y, 'k'), hold on
stem(t, Y2), ylabel('Amplitude'), title('No. of levels, L = 8')
legend('Y_a(t)', 'Y_q[n]')

Y3 = uquant(Y, 16);

43
subplot(2,2,3), plot(t, Y, 'k'), hold on
stem(t, Y3), ylabel('Amplitude'), title('No. of levels, L = 16')
legend('Y_a(t)', 'Y_q[n]')

Y4 = uquant(Y, 32);
subplot(2,2,4), plot(t, Y, 'k'), hold on
stem(t, Y4), ylabel('Amplitude'), title('No. of levels, L = 32')
legend('Y_a(t)', 'Y_q[n]')

The output figure is:

REPORT
[Submit on the next session]
Q1. Write a function quantz such as: [Xq, Eq] = quantz(X, L)where, Xq is the quantized signal, Eq is the
quantization error and L is the number of levels.

Q2. (a)Using the quantz function, quantize the signal: ( ) ( ) for number of levels,
. Take points.
(b) Calculate the mean quantization error, for these levels and plot vs. curve.

44
Reference Books:
(1) “S.K. Mitra”, Digital Signal Processing: A computer based Approach, 3rd Edition.
(2) “J.G. Proakis and D.G. Manolakis”, Digital Signal Processing: Principles, Algorithms and Applications,
4th Edition
Web resource: [Link]

45
 Daffodil International University
Department of Electrical and Electronic Engineering
EEE 322: Digital Signal Processing Laboratory

EXPERIMENT NO: 05

NAME OF THE EXPERIMENT: TIME-DOMAIN ANALYSIS

LTI System Response: Convolution

LTI stands for Linear and Time-Invariant. For a relaxed LTI system, the response ( )to agiven input
signal, ( )can be obtained if we know the system impulse response, ( ). ( )is given by the following relation:

( ) ∑ ( ) ( ) ( ) ( )

which is known as the convolution sum.

Figure 6.1: LTI System Model

Use of conv function

Using the conv function, we can compute the convolution sum of two discrete-time signals. But, the outcome of
conv function does not reveal the timing information. So, time index of the sum should be derived from the signals
to be convolved. Here, the lowest time index of the convolution sum is the sum of the lowest time indices of the two
signals to be convolved. Same idea applies for the highest time index of the sum.

Example 1
Compute the convolution of the two following signals: ( ) 2 3 ( ) 2 3

Solution:
MATLAB Code
x1 = [4 2 6 3 8 1 5];
n1 = -2:4; % generating index
x2 = [3 8 6 9 6 7];
n2 = -4:1;
nmin = n1(1) + n2(1); % left edge of convolved result
nmax = n1(end) + n2(end); % right edge of convolved result
y = conv(x1, x2);
n = nmin:nmax; % generating index of the result
subplot(3,1,1), stem(n1,x1)
subplot(3,1,2), stem(n2,x2)
subplot(3,1,3), stem(n,y)
The result of this example is shown in Fig. 6.2.

46
Exercise 6.1
In example 1, ( ), ( ) and ( ) are all of different lengths and scale and this is far from ideal (Fig. 6.2). Write a
general program to bring the three plots into the same horizontal scale (Fig. 6.3).
m-file which returns both convolved result and index
function [y, ny] = conv_m(x, nx, h, nh)
y = conv(x, h)
nyb = nx(1) + nh(1);
nye = nx(length(x)) + nh(length(h));
ny = linspace(nyb, nye, length(y));

Figure 6.2: With different horizontal scale

Figure 6.3: With same horizontal scale

47
System Response from Difference Equations
An important subclass of discrete-time LTI systems is characterized by the linear constant-coefficient difference
equation:

∑ ( ) ∑ ( ) ( )

where ( ) and ( ) are, respectively, the input and output of the system. * + and * + are the coefficients. We
can modify the equation as:

( ) ∑ ( ) ∑ ( )

( ) ∑ ( ) ∑ ( ) ( )

Use of filter function

Based on ( ), a causal LTI system can be simulated in MATLAB using filter function. Its general form is:
y = filter(b, a, x);
Remember, coefficient must be non-zero. If , filter normalizes the system coefficients by .

Example 2
Consider a system described by the following equation:
( ) ( ) ( )
Find impulse response, step response and sinusoidal response for ( ) ( ) using
filter function. Verify the result using conv function for step and sinusoidal responses.

Solution:
MATLAB Code
a = [1 0.6]; % x(n) coefficient
b = [1 0]; % y(n) coefficient
n = -10:20;

d = [n == 0];
h = filter(b, a, d); % impulse response
subplot(1,3,1), stem(n, h)
title('Impulse response')

u = [n >= 0];
g = filter(b, a, u); % step response
subplot(1,3,2), stem(n, g)
title('Step response')

x = 0.5*sin(n);
y = filter(b, a, x); % sinusoidal response
subplot(1,3,3), stem(n, y)
title('Sinusoidal response')

48
% Verification
y1 = conv(u, h); % step response
y2 = conv(x, h); % sinusoidal response
m = -20:40;
figure(2)
subplot(1,2,1), stem(m, y1)
subplot(1,2,2), stem(m, y2)

Figure 6.4: Impulse, step and sinusoidal response of the system in example 2

Correlation
The cross-correlation function (CCF) of two deterministic real discrete-time finite energy sequences ( ) and ( )
is a third sequence, ( ), defined as:

() ∑ ( ) ( )

In dealing finite duration sequences where ( ) ( ) , ( ) is defined as:

| |

() ∑ ( ) ( )

A distinctpeak in the CCF indicates that the correlated signals are matched for that particular time-shift. Correlation
has important applications in signal detection and system identification.

49
For -point periodic signals, a simpler alternative is:

() ∑ ( ) ( )

Use of xcorr function

Using the xcorr function, we can compute the cross-correlation of two discrete-time signals. But, similar to the
function conv, the xcorr function output does not reveal the timing information. So, time index are derived
from the signals to be correlated. For this reason, the time index of one signal is flipped and negated and then the
sumof the two lowest time indices becomes the lowest time index for the correlated signal. Similarly, the highest
time index can be calculated.

Example 3
Consider, ( ) ( ) where ( ) 2 3. Find and plot the CCF of ( ) and ( ).

Solution:
n = -3:3;
x = [3, 11, 7, 0, -1, 4, 2];
[y, ny] = sigshift(x, n, 2);

rxy = xcorr(x, y);

rxy = rxy/max(rxy);

nyy = -fliplr(ny);
nb = n(1) + nyy(1);
ne = n(end) + nyy(end);
nxy = linspace(nb, ne,length(rxy));

stem(nxy, rxy, 'k')

xlabel('lag, m')
ylabel('r_{xy}(m)')
title('Cross-correlation of x(n) &
y(n)') Figure 6.5: Cross-correlation of ( ) and ( )
From Fig. 6.5, it is apparent that the CCF is very small if , but substantial if . This is because ( )
lags behind ( ) by samples and a match therefore occurs at samples delay.

Using convolution to find CCF: Use of conv

There is a connection between convolution and correlation. For convolution, we set:

( ) ( ) ∑ ( ) ( )

Again, the cross-correlation is defined as:

() ∑ ( ) ( ) ∑ ( ) , ( )- ( ) ( )

So, we can compute CCF using conv function as:

Rxy = conv(x, fliplr(y));
Here, fliplr function folds ( ).

50
Example 3 can also be solved using conv function as follows:

Example 3 (Revisited)
Consider, ( ) ( ) where ( ) 2 3. Find and plot the CCF of ( ) and ( ).
Solution using conv:
MATLAB Code
n = -3:3;
x = [3, 11, 7, 0, -1, 4, 2];

[y, ny] = sigshift(x, n, 2)

[y, ny] = sigfold(y, ny);
[rxy, nxy] = conv_m(x, n, y, ny);
rxy = rxy/max(rxy);

stem(nxy, rxy, 'k')

xlabel('lag, m')
ylabel('r_{xy}(m)')
title('Cross-correlation of x(n) &
y(n)')

Figure 6.6: Estimation of ( )using conv

Note: User-defined functions sigshift, sigfold and conv_m are used in this code. Be sure to save the three functions
in three script files before writing this code.

Auto-correlation
Similar to CCF, the auto-correlation function (ACF) of a signal, ( ) can be defined:

() ∑ ( ) ( ) ( ) ( )

To calculate ACF, just replace ( ) with ( ) in all the previous CCF examples. Note that, ACF of a periodic
signal is also periodic. xcorr function can be used to determine the ACF of a discrete-time signal as:
rxx = xcorr(x);
where, time indices are calculated as before.

Example 4
Find the ACF of the signal:
( ) ( )

Take where period, and ..

Solution: rxx = rxx/max(rxx);

N = 2e-3; nxx = nstep*(1:length(rxx));
nstep = N/100;
n = nstep:nstep:2*N; stem(nxx, rxx, 'k')
x = 2*sin(2*pi*n/N); xlabel('lag, m')
ylabel('r_{xx}(m)')
rxx = xcorr(x); title('Auto-correlation of x(n)')

51
 Figure 6.7: Auto-correlation sequence of ( )

REPORT
[Submit on the next session]
Q1. Consider: ( ) ( ) ( ) ( )
Obtain the step response for this system applying superposition theorem. DO NOT use filter function in that
case. Verify your result using filter function later on.
[Hints: First, find the response, ( )for ( ) ( ), then find the response ( ) for ( ) ( ).
The final response is: ( ) ( ) ( ).]

Q2. Find the ACF of the signal: ( ) . /. Take where period, and
. Explain why ACF magnitude decreases for successive periods?

Q3. Observe the ACF of a random white noise sequence. [Hints: Use the randn function to generate a random
white noise sequence.]

Reference Books:
(1) “S.K. Mitra”, Digital Signal Processing: A computer based Approach, 3rd Edition.
(2) “J.G. Proakis and D.G. Manolakis”, Digital Signal Processing: Principles, Algorithms and Applications,
4th Edition
Web resource: [Link]

52
 Daffodil International University
Department of Electrical and Electronic Engineering
EEE 322: Digital Signal Processing Laboratory

EXPERIMENT NO: 06

NAME OF THE EXPERIMENT: Z-TRANSFORM AND ITS APPLICATION

z-Transform
The z-transform of a discrete time signal is defined as:

( ) ∑ ( )

where, is a complex variable. Since z-transform is an infinite series, it exists only for those values of for which
this series converges. The region of convergence (ROC) of ( ) is defined as the set ofall values of for which
( ) attains a finite value.

Example 1
Let‟s take a finite duration signal like: ( ) 2 3

The z-transformed signal of ( ) is:

( ) ∑ ( )

Its ROC is entire z-plane except .

Again, for ( ) 2 3, the z-transformed signal is:

( ) ∑ ( )

Its ROC is entire z-plane except and .

z-Transform of Different Functions

Region of Convergence
Signal, ( ) z-transformed signal, ( )
(ROC)
( ) 1 All

( ) | |

( ) | | | |

( ) | | | |

( ) | | | |
( )

53
 ( ) | | | |
( )

( ) ( ) | |

( ) ( ) | |

( ) ( ) | | | |

( ) ( ) | | | |

MATLAB Representation
For representing a signal like this ( ) in MATLAB, we just have to put the
coefficients of the polynomial like:
>> X = [1, 2, 5, 7, 0, 1];
For converting the time-domain signal into equivalent -domain, ztrans function can be used. It evaluates signals
of the form ( ) ( ), i.e. for non-negative values of . There is also the iztrans function which performs the
inverse z-transform.

Example 2
To find the z-transform of a function like ( ), write the code:
>> syms a n f;
>> f = a^n;
>> ztrans(f)
Properties of z-transform

54
Exercise 7.1
Let, ( ) and ( ) . Determine ( ) ( ) ( ).
Use conv()function to obtain the result and verify theoretically.

Inverse z-Transform
The inverse z-transform is the conversion of -domain signal intotime-domain signal. The inversion can be done in
many ways such as – Cauchy‟s Integral theorem, long division process, partial fraction expansion and so on.

1. Partial Fraction Expansion

Suppose a function ( ) is given in z-domain( ):
( )
( )
( ) ∑
( )
Performing partial fraction expansion:

( ) ∑[ ] ∑

Then, the time domain signal is:

( ) ∑[ { }] ∑ ( )

MATLAB Representation
The function residuez is available to compute the residue part and the direct terms of a rational part in
( ). The function finds the residues, poles and direct terms of ( ).
[R, p, C] = residuez(b, a)
The two polynomials ( ) and ( ) are given in two vectors b and a, respectively.

Exercise 7.2
A function ( ) is given below:
( )
Find the residue, poles and direct terms of it from the function and hence determine the inverse of ( )
manually. Verify your result using iztrans.

2. Long Division
For causal sequences, the z-transform ( ) can be expanded into power series in . In series expansion, the
coefficient multiplying by is then the nth sample ( ). For a rational ( ), a convenient way to determine
the power series is to express the numerator and denominator as polynomials in and then obtain the power
series expansion by long division.

MATLAB Representation
The function impz is available to provide the time-domain sampled values from ( ).
[h, t] = impz(b, a)
The same process can be performed by using the filter function, with the impulse function as input.
h = filter(b, a, x) % x = delta function

Exercise 7.3
The transfer function of a causal system is given by:
( )
Determine the first 10 coefficients of the impulse response of this system.

55
Pole – zero plot & ROC for different cases
For a system transfer function, ( ), poles are the values for which ( ) becomes infinite. These are the roots of the
denominator of a transfer function, i.e., ( ) and zeros are the values for which ( ) becomes zero. These are the
roots of the numerator of a transfer function, i.e., ( ).
The pole-zero plot of a system provides information about the system behavior.

MATLAB Representation
The function zplane is available in MATLAB for obtaining the pole-zero plot.

Example 3
For a given transfer function:

( )

The pole-zero plot can be obtained

through the code:
>> b = [1 0];
>> a = [1 -0.9];
>> zplane(b, a)

Figure 7.1: Pole-zero plot of a typical transfer function

Region of Convergence (ROC)
For a transfer function of the form:
( )
 The ROC of a causal signal is the exterior of a circle of radius | |.
 The ROC of an anti-causal signal is the interior of a circle of radius | |.
 The ROC of a non-causal signal is a ring (annular region) in the z-plane.

Pole Location & Time-Domain Behavior

Causal real signals with real poles or complex conjugate pair of poles, which are inside or on the unit circle, are
always bounded in amplitude. On the contrary, signals become unbounded with the poles outside of unit circle.

(i) Effect of a real single pole

For a real single pole, the z-transformed function is:
( ) ( ) ( )
The effect of the pole is illustrated in the figure below:

56
 Figure 7.2: Effect of a single real pole position on system response

(ii) Effect of a complex conjugate pair of poles

For a complex conjugate pair of poles, the effect on system response is illustrated below:

(a) Poles: where | | –

(b) Poles: where | | –

57
 (c) Poles: where | | –

(iii) Effect of real multiple poles

For real multiple poles, the z-transformed function is:
( )
( )
The effect is illustrated below:

MATLAB Representation
From the above diagrams, it is evident that different systems provide different impulse responses depending on their
pole-zero location. Hence, different signals can be generated from those defined systems.
Example: A ramp signal can be generated from a system with two positive real poles on the unit circle. Thus, the
following code generates a ramp signal:
>> num = [1 0 0];
>> den = [1 -2 1];
>> [h, yt] = impz(num,den,100);
>> stem(h)

Exercise 7.4
Design a simple low-pass and a high-pass filter.

58
Causality & Stability
An LTI system is –
Causal—If and only if the system transfer function ROC is the exterior of a circle of radius, | | .
BIBO Stable —If and only if the system transfer function ROC includes the unit circle, | | .
So, a causal LTI system is BIBO stable if and only if all the poles of ( )are inside the unitcircle.

MATLAB Representation
The functionzp2sos() can be used to convert the z-domain transfer function into the factored form.
[sos, g] = zp2sos(z, p, k)
Again, thetf2zp() can be used to obtain poles, zeros and gain constant.
[z, p, k] = tf2zp(b, a)
Here, the numerator and denominator polynomials are represented as vectors b and a. Zeros, poles and gain constant
are expressed in z, p and k respectively.

Exercise 7.5
Determine the system function and the response for unit step and unit impulse input described by the difference
equation –
( ) ( ) ( )
(Assume a causal LTI system.)
Looking only the pole-zero plot of ( ), state whether the system will be stable. Comment on ROC.

Exercise 7.6
An LTI system is characterized by the system transfer function –

( )
(a) Express ( ) in factored form and find its poles and zeros.
(b) Specify the ROC of ( ) and determine ( ) for the following conditions:
(i) The system is causal
(ii) The system is anti-causal
(iii) The system is non-causal
For which case, the system is stable?

Pole-Zero Cancellation
When a z-transform has a pole that is at the same location as a zero, the pole is cancelled by zero. Consequently, the
term containing that pole in the term vanishes. This pole-zero cancellation can occur either in the system
function itself or in the product of the system transfer function with the z-transform of the input signal.

Exercise 7.7
Determine the system transfer function and the output response for the causal LTI system described by the
difference equation –

( ) ( ) ( ) ( )

where, ( ) ( ) ( )
Is it possible to reconstruct the original system response from the output response?

59
REPORT
[Submit on the next session]
Q1. Solve exercises 7.1 - 7.7 and answer the questions asked.

Reference Books:
(1) “S.K. Mitra”, Digital Signal Processing: A computer based Approach, 3rd Edition.
(2) “J.G. Proakis and D.G. Manolakis”, Digital Signal Processing: Principles, Algorithms and Applications,
4th Edition.
(3) “V.K. Ingle &J.G. Proakis”, Digital Signal Processing using MATLAB, Edition 2000,
Thomson-Brooks/Cole Co Ltd.
Web resource: [Link]

60
 Daffodil International University
Department of Electrical and Electronic Engineering
EEE 322: Digital Signal Processing Laboratory

EXPERIMENT NO: 07

NAME OF THE EXPERIMENT: FREQUENCY DOMAIN ANALYSIS

Discrete-Time Fourier Series (DTFS)

Consider a periodic sequence ( ) with period , i.e., ( ) ( ) for all . The Fourier series representation
of ( ) consists of harmonically related exponential functions and expressed as:

( ) ∑

Fourier coefficients * + provide the description of ( ) in frequency domain. * + can be

computed as:

∑ ( )

Note that, . That is, * + is a periodic sequence with fundamental period . Thus, the spectrum of a signal
( ), which is periodic with period with period , is again a periodic sequence with period .
Average power can be given as:

∑| | ∑| ( )|

The sequence | | for is the distribution of power as a function of frequency and is called the
Power Density Spectrum (PSD) of the periodic signal.

Calculation of DTFS Coefficients

Example 1
Consider a continuous rectangular pulse sequence of period and a pulse-width of . The signal is
sampled at and sampled discrete-time signal with its index is shown in figure 8.1.

Fig. 8.1: A rectangular pulse in continuous domain and its sampled version

61
Find the DTFS coefficients of the signal and reconstruct the signal from the DTFS coefficients.

Solution:
The MATLAB code for calculating DTFS coefficients is given below:
% Rectangular pulse train...
Fs = 100e3;
dt = 1/Fs; % generating the rectangular pulse train
T = 1e-3; % period of the pulse train
D = 0.1; % duty cycle
PW = D*T; % pulse width
f = 1/T; % analog frequency
t = -T/2:dt:T/2; % time interval for a period
n = t/dt; % index for data points in a period
L = PW/dt; % data points in the the high time
x = zeros(1, length(t)); % initializing a single rectangular pulse
x(find(abs(n) <= L/2)) = 1.1; % generation of single rectangular pulse
figure(1)
subplot(211), plot(t, x, 'k')
title('(a) Original continuous signal')
xlabel('Time (seconds)'), ylabel('x(t)')
subplot(212), stem(n, x, '.k')
title('(b) Sampled discrete signal')
xlabel('Index, n'), ylabel('x(n)')

% DTFS coefficients...
N = length(x); % total no data points in a period
Nc = N; % total no of coefficients
if mod(Nc,2)==0,
k = -Nc/2:Nc/2-1;
else
k = -(Nc-1)/2:(Nc-1)/2;
end
c = zeros(1, length(k)); % initializing fourier coefficients
for i1 = 1:length(k)
for i2 = 1:length(x)
c(i1) = c(i1) + 1/N*x(i2)*exp(-i*2*pi*k(i1)*n(i2)/N);
end
end
figure(2)
subplot(211), stem(k, abs(c), '.k')
title('Fourier series coefficients c_k')
xlabel('k'), ylabel('|c_k|')
subplot(212), stem(k, angle(c)*180/pi, '.k')
xlabel('k'), ylabel('angle(c_k)')

62
figure(3), stem(k*f, c, '.k')
title('Fourier series coefficients c_k')
xlabel('Frequency (Hz)'), ylabel('c_k')

% Reconstruction of signal...
t_remax = T/2;
t_re = -t_remax:dt:t_remax;
n_re = t_re/dt;
x_re = zeros(1, length(n_re));
for i1 = 1:length(k)
for i2 = 1:length(x_re)
x_re(i2) = x_re(i2) + c(i1)*exp(i*2*pi*k(i1)*n_re(i2)/N);
end
end

figure(4)
subplot(211), stem(n_re, x_re, '.k')
title('Reconstructed signal')
xlabel('n'), ylabel('x_{reconstructed}')
subplot(212), plot(t_re, x_re, 'k' ,0, 0)
title('Reconstructed signal')
xlabel('t'), ylabel('x_{reconstructed}')

The output of the code is shown below:

63
Exercise 8.1
Note that in the Example 1, the time interval is defined from –T/2 to T/2,
(i) Change the interval to 0 to T and to –T/4 to 3T/4. Do you observe any change in the magnitude spectrum?
Explain the observation.
(ii) Change the interval to –T to T. Explain the observation.

64
Discrete-Time Fourier Transform (DTFT)
If ( ) is absolutely summable that is

∑ | ( )|

then its discrete time Fourier transform is given by:

( ) ∑ ( )

Now is a real variable between to but one important property of ( ) is its periodicity in with period
. So we only need one period of ( ) i.e., , - , -. The inverse DTFT equation, in other words
the equation for reconstructing the signal from ( ) is given by:

( ) ∫ ( )

Calculation of DTFT Coefficients

Since the frequency is continuous variable, DTFT cannot be implemented in digital hardward in rigorous sense. By
defining a fine grid of frequency in the range , -, the frequency variable is usually handled in DTFT
implementation. In MATLAB, we evaluate ( ) at equidistant frequencies between , - and then interpolate
using plot function.
Say we want equidistant samples between , -. Then –

Then, the DTFT equation is:

( ) ∑ ( )

Example 2
Find the DTFT coefficients of a sinc pulse.

Solution:
The MATLAB code for calculating DTFT coefficients is given below:
% Generating a sinc pulse...
f_c = 1/8; % defining frequency variable for sinc pulse
n = -40:40; % defining the index for sinc pulse
x = sinc(2*f_c*n);
figure(1), stem(n, x, '.k');
title('Discrete time signal')
xlabel('n'), ylabel('x(n)')

% DTFT coefficients...
M = 101; % number of points in digital frequency grid
w = linspace(-pi, pi, M); % defining the digital frequency grid
dw = w(2) - w(1); % resolution of digital frequency
X = zeros(1, M); % initializing the DTFT of x(n)

65
for i1 = 1:M
for i2 = 1:length(x)
X(i1) = X(i1) + x(i2)*exp(-i*w(i1)*n(i2));
end
end

figure(2), plot(w, abs(X), 'k')

xlabel('Frequency (rad/sec)'), ylabel('X(w)')
title('FREQUENCY SPECTRA')

% Reconstruction of signal...
n_re = -80:80;
x_re = zeros(1, length(n_re)); % INITIALIZING RECONSTRUCTED SIGNAL
for i1 = 1:M
for i2 = 1:length(x_re)
x_re(i2) = x_re(i2) + 1/(2*pi)*X(i1)*exp(-i*w(i1)*n_re(i2))*dw;
end
end

figure(3), stem(n_re, x_re, '.k')

xlabel('t'), ylabel('x_{reconstructed}')
title('Reconstructed signal')

The output of the code is shown below:

66
Exercise 8.2
Increase the range of the frequency grid to [0, 2π]. Observe and explain the effects.

Discrete Fourier Transform (DFT)

We know that aperiodic finite energy signals have continuous spectra (DTFT). In case of a finite length sequence
( ) ., Then only values of ( ) over its period, called the frequency samples, are sufficient to
determine ( ) and hence ( ). This leads to the concept of discrete Fourier transform (DFT) which is obtained by
periodic sampling of ( ) (DTFT).
We often compute a higher point ( point) DFT where . This is because padding the sequence ( ) with
zeros and computing an point DFT results in a better display of the Fourier transform ( ). To
summarize, the formulas are (for causal sequence)

( ) ∑ ( ) ∑ ( )

( ) ∑ ( )

In terms of general DFT equation –

Where
, ( ) ( ) ( ) ( )-
, ( ) ( ) ( ) ( )-

( )

( ) ( )
[ ]

Note: matrix can be generated by dftmtx function.

67
Calculation of DFT Coefficients
Since the frequency is continuous variable, DTFT cannot be implemented in digital hardward in rigorous sense. By
defining a fine grid of frequency in the range , -, the frequency variable is usually handled in DTFT
implementation. In MATLAB, we evaluate ( ) at equidistant frequencies between , - and then interpolate
using plot function.
Say we want equidistant samples between , -. Then –

Then, the DTFT equation is:

( ) ∑ ( )

Example 3
Calculate the -point DFT of the sequence, ( ) , -

Solution:
The MATLAB code for calculating DFT coefficients is given below:
% calculation of 5 point DFT
x = [1,1,0, 0, 1];
N = 5; % no of frequency samples
L = length(x);
x = [x zeros(1, N-L)]; % zero padding
n = length(x);
n = 0:N-1; % index of data sequence
k = 0:N-1; % index of frequency sample
Wn = exp(-j*2*pi/N);
WN = Wn.^(n'*k);
X = WN*x'; % DFT

figure(1)
subplot(211), stem(k/N, abs(X), 'k')
title('Magnitude Spectrum')
xlabel('Digital frequency, f (Hz)'), ylabel('|X(f)|')

subplot(212), stem(k/N, angle(X)*180/pi, 'k')

title('Phase Spectrum')
xlabel('Digital frequency, f (Hz)'), ylabel('angle(X(f))')

abs(X)
angle(X)*180/pi

The output of the code is shown below:

>> abs(X) =
3.0000 1.6180 0.6180 0.6180 1.6180
>> angle(X)*180/pi =
0 0 -180.0000 -180.0000 0.0000

68
REPORT
[Submit on the next session]
Q1. Solve the exercises given and answer the questions asked.

Reference Books:
(1) “S.K. Mitra”, Digital Signal Processing: A computer based Approach, 3rd Edition.
(2) “J.G. Proakis and D.G. Manolakis”, Digital Signal Processing: Principles, Algorithms and Applications,
4th Edition.
(3) “V.K. Ingle & J.G. Proakis”, Digital Signal Processing using MATLAB, Edition 2000,
Thomson-Brooks/Cole Co Ltd.
Web resource: [Link]

69