Laboratory Manual

ECE 2108: Signal & Systems Laboratory

1st Edition: 07 May, 2023

Department of Electronics and Communication Engineering

Khulna University of Engineering & Technology
Table of Contents

Experiment 01: Introduction to MATLAB and different singularity functions

generation and plotting. .............................................................................................1

Experiment 02: Study on MATLAB function, shifting of signal, verification of

Euler’s equation, some statistical MATLAB function and solution of linear set of

equation. ...................................................................................................................20

Experiment 03: Study on Fourier series & Fourier transform, correlation &

convolution of signals, solution of differential equations, and parallel & cascaded

systems. ....................................................................................................................30

Experiment 04: Experiment on linearity and time variance of a system, Laplace &

Z transform and stability of a system. ......................................................................55

Experiment 05: Study on Continuous and Discrete State Space Model..................73

Experiment 01: Introduction to MATLAB and different singularity functions
generation and plotting.
1.1. Objectives:
1. To get introduced with MATLAB.
2. To know about different functionality of MATLAB.
3. To plot the signals of some singularity functions on MATLAB such as unit impulse
signal, unit step signal, ramp signal.
4. To plot exponential signal, sinusoidal signals and random signals
1.2. Theory:

MATrix LABoratory is the meaning behind the term MATLAB which was written originally to
provide easy access to matrix software developed by the LINPACK (linear system package) and
EISPACK (Eigen system package) projects. It supports object-oriented programming, has
advanced data structures, and comes with built-in editing and debugging tools. These things make
MATLAB a superior research and teaching tool. Visualization is easy with MATLAB by the use
of some graphics command. There are toolboxes for signal processing, control theory, simulation,
optimization, and several other fields of applied science and engineering.

After opening MATLAB software, the following interface will be seen.

Fig 1.1: User Interface of MATLAB.

1
Toolstrip: The Toolstrip divides the functionality of MATLAB into a number of tabs. Sections of
tabs contain several controls. The controls in MATLAB are the buttons, drop-down menus, and
other user interface components. The Home tab, for example, has sections for operations on files,
variables, code, and other things. There are several controls in the File area, such as writing scripts
(New Script), opening files (Open), and comparing two files (Compare).

Directory of Current Folder: It shows the path of current folder.

Current Folder: For running codes written on each script, the MATLAB file must be kept on
current folder or the path must be added to MATLAB path. Current Folder Window shows the
contents kept on that folder.

Script File: A script is a file that contains multiple sequential lines of MATLAB commands and
function calls. One can run a script by typing its name at the command line or clicking on run icon
from toolstrip. Multiple script file can be opened simultaneously. Script files have a filename
extension .m and are often called M-files.

Command Window: Individual statements can be inserted in the command window at the
command line, indicated by the prompt (>>). Then by pressing “Enter” key it executes the
command and shows the result. The difference between script file and command window is that
script file can contain multiple lines of statement, executes them one by one sequentially and
shows result finally; but command window executes each statement as soon as they are inserted in
the command window and shows the result.

Workspace Window: It shows all the declared variables and their values or dimensions which
MATLAB finds while running a code.

Now, we will see some functionality of MATLAB,

1. Declaring a MATLAB variable: MATLAB variables are created with an assignment

statement and type of the variable does not need to be defined explicitly. Its syntax is,

Variable name = a value (or an expression)

Where expression is a combination of numerical values, mathematical operators, variables, or

function calls. Once a variable is created and the statement has been run, it will be seen the
workspace with its value.

2
2. Suppressing result from Command window: Result of a statement is not shown in the
command window when a semicolon is put at the end of the statement.
3. Controlling the appearance of floating point number: MATLAB by default displays only
4 decimals in the result of the calculations. However, MATLAB does numerical calculations
in double precision, which are 15 digits. The command “format” controls how the results of
computations are displayed.

Fig 1.2: Demonstration of the control of appearance of floating point.

Here, variable x was declared and assigned a value of 759.5791596. But MATLAB showed
only 4 digits after the decimal point. Then using the command “format long”, we can see all
15 digits after the decimal point. To return to the standard format, “format short”, or simply
“format” command is to be used.

4. Clearing workspace: To delete all the stored variables and their values “clear” command is
used.
5. Clearing Command Window: To clear the command window “clc” command is used.
6. Creating matrix: Matrices are the basic elements of the MATLAB environment. A matrix is
a two-dimensional array consisting of m rows and n columns. Special cases are column vectors
(n = 1) and row vectors (m = 1). Here, three rules are to be remembered,
 Matrices are generated with the help of third parenthesis, [ ].
 Columns are separated by the help of comma (,) or space.
 Rows are separated by the help of semicolon (;).

3
If we want to generate a row vector, meaning a matrix with only one row but multiple columns,
all elements of the matrix will be put inside a [ ] parenthesis and the elements will be separated
from each other by the use of comma (,) or space. Comma indicates starting of a new column.
It is shown in Fig 1.3.

Fig 1.3: Generation of row vector.

Now, if we want to generate a column vector, meaning a matrix with only one column but
multiple rows, the elements will be separated from each other by the use of a semicolon (;).
Semicolon indicates the starting of a new row. It is shown in Fig 1.4.

Fig 1.4: Generation of column vector.

Now to create a matrix of m rows and n columns, write all the ‘n’ elements of a row separated
by a comma or space, a new row is started by using semicolon and the process goes on ‘m’
rows and put them all inside a [ ] parenthesis. It is shown in Fig 1.5.

Fig 1.5: Generation of a (2x4) matrix

4
 There are some special built-in functions to generate matrices,

Table 1.1: Special functions to generate matrices

zeros(m,n) Generates a matrix of (m × n) size with all the

elements are 0
ones( m,n) Generates a matrix of (m × n) size with all the
elements are 1
rand(m,n) Generates a matrix of uniformly distributed random
numbers of (m × n) size
randn(m,n) Generates a matrix of random numbers with Gaussian
distribution of (m × n) size
randi([imin, imax], [m,n]) Generates a matrix of integer numbers in between
(imin, imax) of (m × n) size

7. Matrix operation: Some operations of matrices are,

 Transpose of a matrix: Transpose matrix is found by (‘) operator.

Fig 1.5: Finding transpose matrix.

 Size of a matrix: Size of a matrix can be found by the following command,

size(variable)

5
 Accessing Elements: A specific element of a matrix can be found from its location or
indexes. Its command would be,

variable(row_index, column_index)

Indexing for rows and columns starts with 1. So, if we want to get the element of third row
and second column then the command would be,

variable(3,2)

 Accessing Block of Elements: It is possible to get multiple elements from specific rows
and columns. The command is,

Variable ( [indexes of the rows], [indexes of the columns] )

For example, the following command will access the elements of rows 2, 3, 5 and columns
1, 6.

variable ( [2 3 5], [1 6] )

Colon operator (:) can be used for this purpose. It is described later in details. Colon
operator is used to create a row vector within a starting and ending element with a regular
gap. So, if we want to get the elements of rows 1 to 5 and columns of 4,6,8 then the
command will be,

variable ( 1:5, 4:2:8 )

These two commands are shown in Fig 1.6.

6
 Fig 1.6: Finding blocks of elements of a matrix.

Only colon operator with no bound specifies all indexes of row or columns. For example,
the following command will extract elements from all rows but of columns 5 to last.

variable(: , 5:end)

 Inverse of a matrix: It is found by the command,

inv(variable)

 Multiplication of matrix:
- Command for element wise multiplication of A and B matrix: A.*B
- Command for dot multiplication of A and B matrix: A*B

7
  Division of matrix:
- Command for right division of matrix A with matrix B(𝐴𝐵 −1): A/B
- Command for left division of matrix A with matrix B(𝐴−1 𝐵): A\B
- Command for element wise division of matrix A with matrix B: A.\B
8. Colon (:) operator: Colon operator creates a row vector within some specified value with a
certain gap. If the gap is not defined then the gap is by default 1. Its syntax is,

Variable = starting value : gap : ending value

Creating a row vector from 1 to 10 with gap 1 and 2 is shown in Fig 1.7

Fig 1.7: Creating row vector with colon operator.

9. Linear spacing: It also creates a row vector like colon operator but here number of elements
is defined instead of gap. With the provided elements number, gap is automatically calculated
and elements are found. Its syntax is,

variable= linspace(starting value, ending value, number of elements)

8
 Fig 1.8 shows generation of two row vectors of values within 1 to 5 where number of elements
are 3 and 10

Fig 1.8: Creating row vector with linspace operator.

10. Mathematical Functions: MATLAB has many built-in functions to perform mathematical
operation. Some of these functions are described in Table 1.2.

Table 1.2: Some built-in mathematical functions of MATLAB and their description.

Syntax of the function Description

cos(x) Finds cosine of x
sin(x) Finds sine of x
tan(x) Finds tangent of x
acos(x) Finds inverse cosine of x
asin(x) Finds inverse sine of x
atan(x) Finds inverse tangent of x
sqrt(x) Finds square root of x
log(x) Finds natural logarithm of x
log10(x) Finds common logarithm of x
exp(x) Finds exponential of x

9
11. Control Flow: There are 4 control flow operation in MATLAB of which three is described
here,
 if-else: The if...else statement executes two different codes depending upon whether
the test expression is true or false. Its structure is,

if expression
code 1
else
code 2
end

An example code of if-else condition is,

x=10;
if x<5
disp('x is smaller than 5');
else
disp('x is greater than 5');
end

Its output is

Fig 1.9: Output of above if-else code

 for loop: A loop is used to repeat a block of code until the specified condition is met.
Its structure is,

for variable = initial_value:increment: final_value

code 1
end

10
 An example code is,

for i=1:2:6
disp(i)
end

Its output is,

Fig 1.10: Output of above for loop

 while loop: A while loop repeatedly executes a code as long as a given condition is
true. Its structure is,

while expression

code 1

end

An example code can be,

i=1;
while i<6
disp(i)
i=i+2;
end
Its output is,

Fig 1.11: Example of while loop

11
12. Plotting in MATLAB: MATLAB graphing procedure is to take a vector of x-coordinates and
a vector of y-coordinates, locate the points (𝑥𝑖 ,𝑦𝑖 ), with i = 1,2,3…n and then join them by
straight lines. Both x and y coordinate vectors must be of same length. Then the command for
plotting is,

plot(x,y)

And to plot in discrete time domain, following command is used,

stem(x,y)

MATLAB has several built-in functions to make each plot more easy to analyze or visualize.

 figure(number): It finds a figure in which the Number property is equal to n. Then this
figure is made current figure. If no figure exists with that property value, a new figure
is created and Number property is set to n.
 subplot(rcn): In each window or figure there can be multiple subplot which is defined
by a grid of ‘r’ row, ‘c’ column. Then it plots the graph on the position specified by n.
The first subplot is the first column of the first row; the second subplot is the second
column of the first row, and so on. If axes already exist in the specified position, then
this command makes the axes the current axes.
 title(‘txt’): It adds the specified title to the current axis.
 xlabel(‘txt’), ylabel(‘txt’): xlabel(‘txt’) adds labels in the x-axis of the current axes and
ylabel(‘txt’) adds labels in the y-axis.
 hold on: hold on retains plots in the current axes so that new plots added to the axes do
not delete existing plots. New plots use the next colors and line styles based on the
ColorOrder and LineStyleOrder properties of the axes.

Let’s consider we have three random signals of which we want to plot two signals in one
subplot using holdon and the other signal in the remaining subplot. Then the code can be,

clc; clear;
y1=rand(1,10);
y2=rand(1,10);
y3=rand(1,10);
x=1:10;
figure(1)
subplot(121)
plot(x,y1)

12
 hold on
plot(x,y2)
xlabel('x axis');
ylabel('y axis');
title('Plotting example1')
hold off
subplot(122)
plot(x,y3)
xlabel('x axis');
ylabel('y axis');
title('Plotting example2')

Its output will be the Fig 1.12.

Fig 1.12: Plotting of three random signals in two subplot with necessary labeling.

Now with the above information we will try to plot some singularity function in MATLAB.

1. Unit Impulse Function: An unit impulse function is defined as

𝛿(𝑡) = 0, 𝑡≠0
and
+∞
∫ 𝛿(𝑡) 𝑑𝑡 = 1
−∞

These equations state that area of the impulse function is unity and it has approximately
zero width and infinity height. In discrete time domain, this function can be,

1, 𝑛=0
𝛿(𝑛) = {
0, 𝑛≠0

13
Code to visualize Unit Impulse signal:

clc
clear
t = -3:3;
amp = [0,0,0,1,0,0,0];
figure(1); % number of figure
subplot(121)
p=plot(t,amp); % continuous plot
set(p,'Linewidth',2);
title('Unit Impulse Signal(Continuous Time Domain)');
xlabel('Time'); %label of x axis
ylabel('Amplitude'); %label of y axis
subplot(122)
s=stem(t,amp); % discrete plot
set(s,'Linewidth',2);
xlabel('Time'); %label of x axis
ylabel('Amplitude'); %label of y axis
title('Unit Impulse Signal(Discrete Time Domain)');

Output:

Fig 1.13: Plot of Unit Impulse signal

14
2. Unit Step Function: Integration of impulse function generates unit step function.

𝑡
1, 𝑡>0
𝑢(𝑡) = ∫ 𝛿(𝑡) 𝑑𝑡 = {
−∞
0, 𝑡<0

For discrete-time signal,

1, 𝑛≥0
𝑢(𝑛) = {
0, 𝑛<0

Code to visualize Unit Step signal:

clc
clear
t=-5:1:25;
amp=[zeros(1,5), ones(1,26)];
figure(1); % number of figure
subplot(121)
p=plot(t,amp); % continuous plot
set(p,'Linewidth',2);
axis([-5 25 -0.5 1.5]);
title('Unit Step Signal(Continuous Time Domain)');
xlabel('Time'); %label of x axis
ylabel('Amplitude'); %label of y axis
subplot(122)
s=stem(t,amp); % discrete plot
axis([-5 25 -0.5 1.5]);
set(s,'Linewidth',2);
xlabel('Time'); %label of x axis
ylabel('Amplitude'); %label of y axis
title('Unit Step Signal(Discrete Time Domain)');

Output:

Fig 1.14: Plot of Unit Step signal

15
3. Unit Ramp Function: Double integration of unit impulse function or single integration
of unit step signal generates unit ramp function.

𝑡 𝛼
𝑟(𝑡) = ∫ ∫ 𝛿(𝑡)𝑑𝑡 𝑑𝛼
−∞ −∞

𝑡
𝑜𝑟, 𝑟(𝑡) = ∫ 𝑢(𝛼)𝑑𝛼
−∞

𝑡
𝑜𝑟, 𝑟(𝑡) = ∫ 1𝑑𝛼
0

𝑡, 𝑡>0
𝑜𝑟, 𝑟(𝑡) = {
0, 𝑡<0

For discrete time domain,

𝑛, 𝑛≥0
𝑟(𝑛) = {
0, 𝑡<0

Code to visualize Unit Ramp signal:

clc
clear
t=-10:10;
amp=[zeros(1,10) t(11:21)];
figure(1); % number of figure
subplot(121)
p=plot(t,amp); % continuous plot
set(p,'Linewidth',2);
axis([-5 10 -0.5 11]);
title('Unit Ramp Signal(Continuous Time Domain)');
xlabel('Time'); %label of x axis
ylabel('Amplitude'); %label of y axis
subplot(122)
s=stem(t,amp); % discrete plot
axis([-5 10 -0.5 11]);
set(s,'Linewidth',2);
xlabel('Time'); %label of x axis
ylabel('Amplitude'); %label of y axis
title('Unit Ramp Signal(Discrete Time Domain)');

16
 Output:

Fig 1.15: Plot of Unit Ramp signal

Now, we will try to plot exponential signal, sinusoidal signal and random signal within one figure,
Code:
clc,
clear

%exponential signal
a= 2; %input('Enter The value of a:');
t=-5:0.5:5;
amp=exp(-a*t)+exp(a*t);
figure(1)
subplot('131');
p=plot(t,amp);
set(p,'Linewidth',2);
title('Exponential Signal');
xlabel('Time');
ylabel('Amplitude');

% Sinusoidal Signal
A= 5; %input('Enter the amplitude:');
f= 1; %input('Enter the frequency:');
t=linspace(1,10,1000);
amp=A*sin(2*pi*f*t);
subplot('132');
pp=plot(t,amp);
set(pp,'Linewidth',2);
title('Sinusoidal Signal');

17
 xlabel('Time');
ylabel('Amplitude');

%random signal
t=-10:1:20;
amp=rand(1,31);
subplot('133');
ppp=plot(t,amp);
set(ppp,'Linewidth',2);
title('Random Signal');
xlabel('Time');
ylabel('Amplitude');

Output:

Fig 1.16: Plot of exponential, sinusoidal and random signal

18
1.3. Experimental Task:
1.

19
Experiment 02: Study on MATLAB function, shifting of signal, verification of
Euler’s equation, some statistical MATLAB function and solution of linear set of
equation.
2.1. Objectives:
1. To know how to define a function in MATLAB.
2. To find even and odd component of a signal using function.
3. To perform shifting operation on some signal.
4. To plot two sinusoidal signal of 90° out of phase.
5. To verify Euler’s equation.
6. To know about some built-in statistical function of MATLAB.
7. To find the solution of a set of linear equation in MATLAB.
2.2. Theory:
Function: Function contains a set of statements which are executed when they are called.
A function can be called multiple times to provide reusability. By using functions, we can
avoid rewriting same code again and again in a program.

In MATLAB, functions are defined in separate files. The name of the file and of the
function should be the same. Its structure is,

function [out1, out2, …, out n]=functionName(ip1,ip2,..,ip M)

statement1;
statement2;
end

Code 1: An example of the definition of function to convert degree into radian can be,

function y = my_function(x)
y=x*(2*pi/360);
end

Write this code in a separate MATLAB file and save it with same name as the function
name. Then open a new MATLAB script file and write a code to call the function. Save
this file in the same folder which contains the file of MATLAB function. It should be
like,
phase=my_function(90);
fprintf('Phase in radian=%f\n', phase)

20
Run this code and in command window following result will be shown,

Fig 2.1: Output of the demo function code to convert degree into radian.

Code 2: A code to find the even and odd components using function of the signal,

𝑥 = 𝑠𝑖𝑛𝑡 + 𝑐𝑜𝑠𝑡

Ultimately, this function is, 𝑥 = √2 sin(𝑡 + 45°) and it’s even component will be 𝑐𝑜𝑠𝑡
and odd component will be, 𝑠𝑖𝑛𝑡. Let’s verify this signals.

Code for function:

function [xe,xo] = components(t,x)
x_reverse = fliplr(x); % time reversal
xe = 0.5*(x + x_reverse); %even component
xo = 0.5*(x - x_reverse); %odd component
subplot(3,1,1);
plot(t,x);
title('Your signal x')
subplot(3,1,2);
plot(t,xe);
title('Even part');
subplot(3,1,3);
plot(t,xo);
title('Odd part');
end

Code for signal generation and calling the function:

t=linspace(0,10,1000);
fs=1;
x=sin(2*pi*fs*t)+cos(2*pi*fs*t);
components(t,x);

21
Its output is,

Fig 2.2: Even and odd components of signal, 𝑥 = 𝑠𝑖𝑛𝑡 + 𝑐𝑜𝑠𝑡

Shifting: A signal x(n) can be either advanced or delayed in time axis. The shifted signal
is represented by, x(n-k).

 If k is positive, the signal is delayed in time axis by ‘k’ unit (Right shift).
 If k is negative, the signal is advanced in time axis by ‘k’ unit (Left Shift).

22
Code 3: Shift a unit step signal by any unit in either right or left direction.

Code for function:

function y = shift(x,k)
if(k>0)
y = [zeros(1,abs(k)) x(1:end-abs(k))];
else
y = [x((abs(k)+1):end) zeros(1,abs(k))];
end
end
Code for signal generation and calling the function:

n=-5:1:25;
x=[zeros(1,5), ones(1,26)];
k=input('Enter Shifting unit:')
y=shift(x,k);
figure(1)
subplot(211)
stem(n,x)
title('Unit step signal')
subplot(212)
stem(n,y)
title('Shifted Unit step signal')

Fig 2.3: Plot of shifted unit step signal by 3 unit

23
Code 4: Plot of two sinusoidal of 90 degree out of phase

clc;
clear all;
fs=1000; % Assume sampling interval
f=2; % frequency of sine wave
t=0:1/fs:1; % Generate Time Vector
phase=90*(2*pi/360); % 90 degree phase convert to radian
x1=sin(2*pi*f*t);
x2=sin(2*pi*f*t-phase);
hold on; %Retain current plot when adding
new plots
plot(t,x1,'k');
plot(t,x2,'r'); grid
xlabel('Tme(sec.)');
ylabel('x(t)');
legend('x1(t)','x2(t)'); %Add legend to graph

Its output is,

Fig 2.4: Plot of sinusoidal signal with 90° out of phase

24
Euler’s equation: Euler's formula relates the complex exponential to the cosine and sine
functions. Its equation is,

𝑒 𝑗𝑡 = 𝑐𝑜𝑠𝑡 + 𝑗 𝑠𝑖𝑛𝑡

Now we plot the real and imaginary part of this complex exponential and verify this
equation.

Code 5: Verification of Euler’s equation:

clc;
clear all;
fs=1000; % Assume sampling interval
t=0:1/fs:1; % Generate Time Vector
z=exp(j*2*pi*t);
plot(t,real(z),'k',t, imag(z), 'r'); grid
xlabel('Time(sec.)');
ylabel('x(t)');
legend('Real x(t)','Imag x(t)');
Its output is,

Fig 2.5: Plot of real and imaginary terms of complex exponential

From the figure it can be seen that the real term of complex exponential is a cosine wave
and imaginary term is a sine wave. Hence, Euler’s formula is verified.

25
Some built in statistical function of MATLAB: There are some built-in statistical
functions in MATLAB, some of them are listed in Table 4.1.

Table 4.1: List and description of some statistical functions.

Functions Description
min(x) Finds minimum value of x
max(x) Finds maximum value of x
mean(x) Finds average value of x
var(x) Finds variance of x
std(x) Finds the standard deviation of x
sqrt(x) Finds the square root of x

Code 6: Code to demonstrate these functions:

clc;
clear all;
fs=1000;
f=2;
t=0:1/fs:1;
x=sin(2*pi*f*t);
xmin=min(x)
xmax=max(x)
xmean=mean(x)
xvar=var(x)
xstad=std(x)
RMS=sqrt(mean(x.^2))

Its output is,

26
Fig 2.6: Operations of these statistical functions on sin(2𝜋 ∗ 2 ∗ 𝑡)

27
Solution of Linear set of equation: Let a set of linear equation is,

𝑎1 𝑥 + 𝑎2 𝑦 + 𝑎3 𝑧 = 𝑝1
𝑏1 𝑥 + 𝑏2 𝑦 + 𝑏3 𝑧 = 𝑝2
𝑐1 𝑥 + 𝑐2 𝑦 + 𝑐3 𝑧 = 𝑝3

Let’s express them is coefficient matrix, variable matrix and constant matrix,
𝑎1 𝑎2 𝑎3
Coefficient matrix, 𝐴 = [ 𝑏1 𝑏2 𝑏3 ]
𝑐1 𝑐2 𝑐3

𝑥
𝑦
Variable matrix, 𝑋 = [ ]
𝑧

𝑝1
𝑝
Constant matrix, 𝐵 = [ 2 ]
𝑝3

Then the solution can be found as,

𝑋 = 𝐴−1 𝐵

Code 7: Code to find solution of linear set of equation:

clc;
clear all;
A=[1 0 0; 2 2 0; 3 3 3];
b=[2 2 2]';
x=inv(A)*b;
disp('Ax=b')
A
b
x

28
 Its output is,

Fig 2.7: Finding solutions of a linear set of equation

2.3. Experimental Tasks:

1.

29
Experiment 03: Study on Fourier series & Fourier transform, correlation &
convolution of signals, solution of differential equations, and parallel & cascaded
systems.
3.1. Objectives:
1. To represent a pulse train with Fourier series and plot on MATLAB and compare the
signal with ideal one.
2. To represent a triangular wave with Fourier series and plot on MATLAB and compare
the signal with ideal one.
3. To represent an exponential signal with Fourier series and plot on MATLAB and
compare the signal with ideal one.
4. To find Fourier transform of single frequency sine wave, double frequency sine wave
and noisy double frequency sine wave.
5. To find Fourier transform of square pulse.
6. To find correlation between sine, cosine and square wave.
7. To find system response to an input 𝑥(𝑡) using convolution operation.
8. To find the solution of a differential equation which represents the system.
9. To find impulse and step response of a system which represented by a differential
equation.
10. To find the output response for a particular input in cascaded system and parallel
connected system.

3.2. Theory:

Any periodic function 𝑓(𝑡) can be expressed in terms of a trigonometric series named
Fourier series,
∞ ∞
𝑎0
𝑓(𝑡) = + ∑ 𝑎𝑛 cos 𝑛𝜔0 𝑡 + ∑ 𝑏𝑛 sin 𝑛𝜔0 𝑡
2
𝑛=1 𝑛=1

Where,
𝑇
2 2
𝑎0 = ∫ 𝑓(𝑡) 𝑑𝑡
𝑇 −𝑇
2

30
 𝑇
2 2
𝑎𝑛 = ∫ 𝑓(𝑡) cos 𝑛𝜔0 𝑡 𝑑𝑡
𝑇 −𝑇
2
𝑇
2 2
𝑏𝑛 = ∫ 𝑓(𝑡) sin 𝑛𝜔0 𝑡 𝑑𝑡
𝑇 −𝑇
2

And, 𝜔0 is fundamental frequency.

3.2.1 Pulse Train: The equation of a pulse train with Fourier series representation is:
𝑛
𝐴 2𝐴 1
𝑓(𝑡) = + ∑ sin(𝑖 ∗ 𝜔0 ∗ 𝑡)
2 𝜋 𝑖
𝑖=1,
𝑖=𝑖+2

Code:

clc; clear all;

N=1000;
t=linspace(0,1,N);
A=5;
res=A/2;
M=input('Number of harmonics=')

for i=1:2:M
e=(2*A/pi)*(1/i)*sin(i*2*pi*1*t);
res=res+e;
subplot(511)
plot(t,e);
hold on
subplot(512)
plot(t,res);
hold on
end
hold off
subplot(511)
ylabel('Amplitude')
xlabel({'time'})
title('sinusoidal signals')
subplot(512)
ylabel('Amplitude')
xlabel({'time'})
title('summation of sinusoidal signals')

subplot(513)
plot(t,res,'k','linewidth',2)
ylabel('Amplitude')
xlabel({'time'})
title(['Resulatant pulse train with ' num2str(M) ' harmonics'])

%ideal pulse train.

31
y1=A*ones(1,N/2);
y2=zeros(1,N/2);
y=[y1 y2];
subplot(514)
plot(t, y, 'k','linewidth',2)
ylabel('Amplitude')
xlabel({'time'})
ylim([-1 A+1]);
title('Ideal pulse train')

%error signal
subplot(515)
plot(t, (res-y), 'k','linewidth',2)
ylabel('Amplitude')
xlabel({'time'})
title('Error signal')

This code will ask for required number of harmonics. For larger harmonics
number, error between ideal and generated signal will reduce.

Fig 3.1: Pulse train generated with 10 harmonics

32
 Fig 3.2: Pulse train generated with 2000 harmonics

It can be seen that with 2000 harmonics the error was very small.

33
3.2.2 Triangular wave: The equation of a pulse train with Fourier series representation
is:
𝑛
8𝐴 1
𝑓(𝑡) = 2 ∑ 2 sin(𝑗 ∗ 𝜔0 ∗ 𝑡)
𝜋 𝑗
𝑗=1
𝑗=𝑗+2

Code:

for i=1:4;
N=1000; % Number of points in one cycle
t=linspace(0,1,N);
sum=0;
M=2^i; % How many harmonics used in particular plot
sign=-1; % 'sign' variables to keep track the sign of sine
term
for j=1:2:M;
sign=sign*(-1);
sum=sum+(sign*8)/(pi*pi)*sin(2*pi*j*t)/(j^2);
end

subplot(2,2,i); plot(t,sum,'r');
title(['Number of Harmonics used:',num2str(M)]);
hold on

%ideal case
y1=linspace(0,1,N/4); % First (1/4)th signal
y2=linspace(1,-1,N/2); % Middle (1/2) signal
y3=linspace(-1,0,N/4); % Last (1/4)th signal
y=[y1 y2 y3];
subplot(2,2,i); plot(t,y,'g','Linewidth',2)

%error signal
subplot(2,2,i);
error=abs(sum-y);
plot(t,error,'k');
axis([0 1 -1 1]);grid
legend('Given shape','Fourier representation','Error');
end

34
 Its output is,

Fig 3.3: Plot of triangular wave and comparison with ideal signal.

3.2.3 Exponential signal: The equation of a exponential signal with Fourier series
representation
𝑛
0.504 ∗ 2
𝑓(𝑡) = 0.504 + ∑ (cos 2𝑗𝑡 + 4𝑗 sin 2𝑗𝑡)
1 + 16𝑗 2
𝑗=1

Code:

plotindex=1;
for i=4:7
N=1000; % Number of points in one cycle
t=linspace(0,pi,N);
sum=0.504;
M=2^i; % How many harmonics used in particular plot
for j=1:M;

sum=sum+((0.504*2)/(1+16*j^2))*(cos(2*j*t)+4*j*sin(2*j*t));
end
subplot(2,2,plotindex); plot(t,sum,'r');
title(['Number of Harmonics used:',num2str(M)]);
hold on

35
%ideal case
y=exp(-.5*t);
subplot(2,2,plotindex);
plot(t,y,'g','Linewidth',2)
hold on

subplot(2,2,plotindex);
error=abs(sum-y);
plot(t,error,'k');
xlim([0,3]);grid
legend('Fourier representation','Ideal shape','Error');
plotindex=plotindex+1;
end

Its output is;

Fig 3.4: Plot of exponential signal and comparison with ideal signal.

36
Fourier transform: Fourier transform converts a time domain signal into frequency
domain i.e. it finds the frequency component of a signal.

3.2.4 Fourier transform of a single frequency sine wave:

Code:

Fs = 1000; % Sampling frequency

T = 1/Fs; % Sample time
t = 0:1/Fs:1;
L=length(t)

% a 50 Hz sinusoid
y = 0.7*sin(2*pi*50*t);
subplot(311)
plot(t(1:50),y(1:50))
title('Single frequency Sine Signal');
xlabel('time (seconds)');

NFFT = 2^nextpow2(L); % Next power of 2 just greater than length

of y
Y = fft(y,NFFT)/L;
Abs_Y = abs(Y);
subplot(312)
plot(Abs_Y);
xlabel('Frequency (Hz)');
ylabel('|Y(f)|');
title('Double-Sided Amplitude Spectrum of y(t)');

% Plot single-sided amplitude spectrum.

f = (Fs/2)*linspace(0,1,NFFT/2);
subplot(313)
plot(f,2*(Abs_Y(1:NFFT/2))) % Multiplication by 2, as amplitude
is splitted into 2 spectrums
title('Single-Sided Amplitude Spectrum of y(t)');
xlabel('Frequency (Hz)');
ylabel('|Y(f)|');

Its output is,

37
 Fig 3.5: Fourier transform of 𝑦 = 0.7 ∗ sin(2𝜋 ∗ 50𝑡).

3.2.5 Fourier transform of sinusoidal signal of two frequency:

Code:

Fs = 1000; % Sampling frequency

T = 1/Fs; % Sample time
t = 0:1/Fs:1;
L=length(t)

y = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
subplot(311)
plot(t(1:50),y(1:50))
title('Double frequency Sine Signal');
xlabel('time (seconds)');

NFFT = 2^nextpow2(L); % Next power of 2 just greater than length

of y
Y = fft(y,NFFT)/L;
Abs_Y = abs(Y);
subplot(312)
plot(Abs_Y);
xlabel('Frequency (Hz)');
ylabel('|Y(f)|');
title('Double-Sided Amplitude Spectrum of y(t)');

% Plot single-sided amplitude spectrum.

f = (Fs/2)*linspace(0,1,NFFT/2);
subplot(313)
plot(f,2*(Abs_Y(1:NFFT/2)))
title('Single-Sided Amplitude Spectrum of y(t)');

38
 xlabel('Frequency (Hz)');
ylabel('|Y(f)|');

Its output is,

Fig 3.6: Fourier transform of 𝑦 = 0.7 ∗ sin(2𝜋 ∗ 50𝑡) + sin(2𝜋 ∗ 120𝑡).

3.2.6 Fourier transform of noisy sinusoidal signal of two frequency:

Code:

Fs = 1000; % Sampling frequency

T = 1/Fs; % Sample time
t = 0:1/Fs:1;
L=length(t)

y = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t)+randn(size(t));
subplot(311)
plot(t(1:50),y(1:50))
title('Noisy Double frequency Sine Signal');
xlabel('time (seconds)');

NFFT = 2^nextpow2(L); % Next power of 2 just greater than length

of y
Y = fft(y,NFFT)/L;
Abs_Y = abs(Y);
subplot(312)
plot(Abs_Y);
xlabel('Frequency (Hz)');
ylabel('|Y(f)|');
title('Double-Sided Amplitude Spectrum of y(t)');

39
% Plot single-sided amplitude spectrum.
f = (Fs/2)*linspace(0,1,NFFT/2);
subplot(313)
plot(f,2*(Abs_Y(1:NFFT/2))) % Multiplication by 2, as amplitude
is splitted into 2 spectrums
title('Single-Sided Amplitude Spectrum of y(t)');
xlabel('Frequency (Hz)');
ylabel('|Y(f)|');

Its output is,

Fig 3.7: Fourier transform of noisy 𝑦 = 0.7 ∗ sin(2𝜋 ∗ 50𝑡) + sin(2𝜋 ∗ 120𝑡).

40
3.2.7 Fourier transform of square pulse:

Code:

Fs = 1000;
T = 1/Fs;
t = 0:1/Fs:1;
L=length(t)

y = [ones(1,10) zeros(1,990)];
subplot(211)
plot(t(1:50),y(1:50))
title('Single frequency Sine Signal');
xlabel('time (seconds)');

NFFT = 2^nextpow2(L);
Y = fft(y,NFFT)/L;
Abs_Y = abs(Y);
Abs_Y = fftshift(Abs_Y);
f = (Fs/2)*linspace(-1,1,NFFT);
subplot(212)
plot(f,2*(Abs_Y(1:NFFT)));
title('Amplitude Spectrum of y(t)');
xlabel('Frequency (Hz)');
ylabel('|Y(f)|');

Its output is,

Fig 3.8: Fourier transform of square pulse.

41
Correlation: Correlation measures the similarity between two data sequences.
Applications: speech processing, image processing, sonar and radar.

3.2.8 The correlation of sin wave and cosine wave and square wave:

Code:

N=1000; % number of points

fs=1000; % sampling frequency
t=(0:N-1)/fs; % time vector
f=1;

x=sin(2*pi*f*t); % 1 Hz sine wave

y=cos(2*pi*f*t); % 1 Hz cosine wave
z=[ones(1,N/2),-ones(1,N/2)]; % 1 Hz square wave
[rxx,lag1]=xcorr(x,x);
[rxy,lag2]=xcorr(x,y);
[rxz,lag3]=xcorr(x,z);

subplot(211)
plot(t,x,'g');
hold on
plot(t,y,'b');
plot(t,z,'k');
legend('sine wave','cosine wave','square wave');
hold off
subplot(212)
plot(lag1,rxx, 'g');
hold on;
plot(lag2,rxy, 'b');
plot(lag3,rxz, 'k'); grid
hold off
legend('rxx','rxy','rxz');

42
 Its output is,

Fig 3.9: Correlation of sine, cosine and square wave.

Convolution: The output of a system for an input expressed as weighted superposition.

Let, a system with impulse response ℎ(𝑡), then, the system can be,

Fig 3.10: System impulse response

For this system with input 𝑥(𝑡), output 𝑦(𝑡) will be,
∞
𝑦(𝑡) = 𝑥(𝑡) ∗ ℎ(𝑡) = ∫ 𝑥(𝜏)ℎ(𝑡 − 𝜏)𝑑𝜏
−∞

3.2.9 We have to find system output if

1, 0≤𝑡≤4
𝑥(𝑡) = {
0, 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒

And

1, 0≤𝑡≤4
ℎ(𝑡) = {
0, 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒

43
Then convolution of the two signal can be found by MATLAB. Its input and output
by mathematical calculation is,

Fig 3.11: Signal of 𝑥(𝑡) and ℎ(𝑡).

Fig 3.12: Signal of 𝑦(𝑡).

Code for function:

function[y,ty]=convint(x,tx,h,th)
%Inputs:
%x is the input signal vector
%tx is the times of the samples in x
%h is the impulse response vector
%th is times of the samples in h
%outputs:
%y is the output signal vector,
%length(y)=length(x)+length(h)-1
%ty is the time of the samples in y
dt=tx(2)-tx(1);
y=conv(x,h)*dt;
ty=(tx(1)+th(1))+[0:(length(y)-1)]*dt;

Code for finding convolution:

clc;
clear all;
close all;
tx=[0:0.01:4];
x=ones(1,length(tx));
th=[0:0.01:4];
h=ones(1,length(th));
[y ty]=convint(x,tx,h,th);
figure;
plot(ty,y); xlabel('Time'); ylabel('Amplitude');axis([0 8 0 4]);

44
 Its output is,

Fig 3.13: Coded convolutional output of the problem.

3.2.10 Solving differential equation: A general Nth-order linear constant-coefficient
differential equation is given by,
𝑁 𝑀
𝑑𝑛 𝑦(𝑡) 𝑑 𝑛 𝑥(𝑡)
∑ 𝑎𝑛 = ∑ 𝑏𝑛 (1)
𝑑𝑡 𝑛 𝑑𝑡 𝑛
𝑛=0 𝑘=0

Then the solution or response 𝑦(𝑡) can be found by the command,

𝑑𝑠𝑜𝑙𝑣𝑒(‘𝑒𝑞𝑛1’, ‘𝑒𝑞𝑛2’, . . . )
Which accepts symbolic equations representing ordinary differential equations and
initial conditions.
For example, we have to find the response of RLC circuit shown in Fig 3.13,

Fig 3.14: Circuit diagram of a series RLC circuit.

45
 Where, R =3 Ohms, L =1 Henry, and C = 12Farad
Its corresponding differentiation equation will be,

𝑑2 𝑦(𝑡) 𝑑𝑦 𝑑𝑥
2
+3 + 2𝑦(𝑡) =
𝑑𝑡 𝑑𝑡 𝑑𝑡

Code for solution:

y = dsolve('D2y+3*Dy+2*y=-30*exp(-3*t)', 'y(0)=0', 'Dy(0)=5',

't');
disp (['y(t) = (', char(y), ')u(t)']);

Its output is,

Fig 3.15: Solution of the differential equation.

3.2.11 Determining Impulse Response and Step Response for a Linear System
Described by a Differential Equation Using MATLAB

For a system described by Eq. (1), the impulse response can be determined by using
the following MATLAB command:

ℎ = 𝑖𝑚𝑝𝑢𝑙𝑠𝑒(𝑏, 𝑎, 𝑡)

For the above described system and differential equation,

𝑑2 𝑦(𝑡) 𝑑𝑦 𝑑𝑥
2
+3 + 2𝑦(𝑡) =
𝑑𝑡 𝑑𝑡 𝑑𝑡

Impulse and step response can be found by coding. However this will require the
above described convolutional function ‘convint’ to be implemented for
finding step response.

46
Code:

clc;
clear all;
close all;
th=0:.01:4;
b = [1];
a = [1 3 2];
h=impulse(b,a,th);
tx=[0:0.01:4];
x=[ones(1,length(tx))];
[y ty]=convint(x,tx,h,th);
subplot(211)
plot(th,h)
xlabel('Time')
ylabel('Amplitude');
title('Impulse response of the system')
subplot(212)
plot(ty,y);
xlabel('Time')
ylabel('Amplitude');
title('Step response of the system')

Its output is,

Fig 3.16: Impulse and step response of the series RLC circuit.

47
3.2.12 System output for Cascaded system

A series cascaded system is shown in figure 3.17.

Fig 3.17: Figure of a cascaded system

For this system, output is

𝑦(𝑡) = [𝑥(𝑡) ∗ ℎ1 (𝑡)] ∗ ℎ2 (𝑡)

Let,

0, 0≤𝑡≤3
𝑥(𝑡) = { 1, 3≤𝑡≤5
0, 5 ≤ 𝑡 ≤ 10

ℎ1 (𝑡) = 3𝑒𝑥𝑝−3𝑡

ℎ2 (𝑡) = 2𝑒𝑥𝑝−2𝑡

Code (convint function is required):

clc;
clear all;
close all;
tx=[0:0.01:10];
tx1=[0:0.01:3];
tx2=[0:0.01:5];
x=[zeros(1,length(tx1)) ones(1,(length(tx)- length(tx1)-length(tx2))) zeros(1,length(tx2))];
th=[0:0.01:10];
h1 =3* exp(-3*th);
[y1 ty1]=convint(x,tx,h1,th);
h2 =(2)* exp(-2*th);
[y ty]=convint(y1,ty1,h2,th);
%using equivalence
[he, the]=convint(h1,th,h2,th);
[yy, tyy]=convint(x,tx,he,the);
%%%%
figure;
subplot(611)
plot(tx,x)
xlabel('Time')
ylabel('Amplitude');

48
title('input x(t)')
ylim([-.05 1.5])
subplot(612)
plot(th,h1)
title('impulse response of first system, h1(t)')
xlabel('Time')
ylabel('Amplitude');
subplot(613)
plot(ty1,y1);
xlabel('Time')
ylabel('Amplitude');
title('intermediate output, y1(t)')
subplot(614)
plot(th, h2)
title('impulse response of second system, h2(t)')
xlabel('Time')
ylabel('Amplitude');
subplot(615)
plot(ty,y)
xlabel('Time')
ylabel('Amplitude');
title('final output, y(t)')
subplot(616)
plot(tyy,yy)
xlabel('Time')
ylabel('Amplitude');
title('final output(using equiavalence), y(t)')

Its output is

49
Fig 3.18: Output response of cascaded system

50
3.2.13 System output for parallel connected system

A parallel connected system is shown in figure 3.19.

Fig 3.19: Figure of a parallel system

For this system, output is

𝑦(𝑡) = 𝑥(𝑡) ∗ [ℎ1 (𝑡) + ℎ2 (𝑡)]

Let,

0, 0≤𝑡≤3
𝑥(𝑡) = { 1, 3≤𝑡≤5
0, 5 ≤ 𝑡 ≤ 10

ℎ1 (𝑡) = 3𝑒𝑥𝑝−3𝑡

ℎ2 (𝑡) = 0.5𝑒𝑥𝑝−0.5𝑡

Code (convint function is required):

clc;
clear all;
close all;
tx=[0:0.01:10];
tx1=[0:0.01:3];
tx2=[0:0.01:5];
x=[zeros(1,length(tx1)) ones(1,(length(tx)- length(tx1)-length(tx2))) zeros(1,length(tx2))];
th=[0:0.01:10];
h1 =3* exp(-3*th);
h2 =0.5* exp(-0.5*th);
[y1 ty1]=convint(x,tx,h1,th);
[y2 ty2]=convint(x,tx,h2,th);
y=y1+y2;
%using equivalence
h=h1+h2;
[yy, tyy]=convint(x,tx,h,th);

51
%%%
figure;
subplot(611)
plot(tx,x)
ylabel('Amplitude');
title('input x(t)')
ylim([-.05 1.5])
subplot(712)
plot(th,h1)
title('impulse response of first system, h1(t)')
ylabel('Amplitude');
subplot(713)
plot(ty1,y1);
ylabel('Amplitude');
title('output of first system, y1(t)')
subplot(714)
plot(th, h2)
title('impulse response of second system, h2(t)')
ylabel('Amplitude');
subplot(715)
plot(ty2,y2)
ylabel('Amplitude');
title('output of second system, y2(t)')
subplot(716)
plot(ty2,y)
ylabel('Amplitude');
title('final output, y(t)')
subplot(717)
plot(tyy,yy)
ylabel('Amplitude');
title('final output(using quivalent circuit), y(t)')

Its output is,

52
Fig 3.20: Output response of parallel connected system

53
3.3. Experimental Task:

54
Experiment 04: Experiment on linearity and time variance of a system, Laplace &
Z transform and stability of a system.
4.1. Objectives:
1. To check whether a system is linear or not.
2. To check whether a system is time-invariant or not.
3. To perform Laplace and inverse Laplace transform.
4. To perform Z and inverse Z transform.
5. To check whether a system is stable or not.
4.2. Theory:

Linear system: A linear system is one in which the principle of superposition and
homogeneity holds.

For a system with two inputs 𝑥1 (𝑡) and 𝑥2 (𝑡) the principle of superposition is:

𝐻[𝑥1 (𝑡) + 𝑥2 (𝑡)] = 𝐻[𝑥1 (𝑡)] + 𝐻[𝑥2 (𝑡)]

Where,

𝐻[𝑥1 (𝑡)] = 𝑜𝑢𝑡𝑝𝑢𝑡 𝑟𝑒𝑠𝑝𝑜𝑛𝑠𝑒 𝑓𝑜𝑟 𝑖𝑛𝑝𝑢𝑡 𝑥1 (𝑡) = ℎ(𝑡) ∗ 𝑥1 (𝑡) = 𝑦1 (𝑡)

The principle of homogeneity is:

𝐻[𝑎1 𝑥1 (𝑡)] = 𝑎1 𝐻[𝑥1 (𝑡)]

𝐻[𝑎2 𝑥2 (𝑡)] = 𝑎2 𝐻[𝑥2 (𝑡)]

Where 𝑎1 & 𝑎2 are weights.

Thus the condition for Linearity will be,

𝐻[𝑎1 𝑥1 (𝑡) + 𝑎2 𝑥2 (𝑡)] = 𝑎1 𝐻[𝑥1 (𝑡)] + 𝑎2 𝐻[𝑥2 (𝑡)]

It is illustrated in Fig 4.1,

55
 Fig 4.1: Illustration of Superposition principle

Thus, a linear system is defined as the system whose,

𝑹𝒆𝒔𝒑𝒐𝒏𝒔𝒆 𝒕𝒐 𝒕𝒉𝒆 𝒔𝒖𝒎 𝒐𝒇 𝒘𝒆𝒊𝒈𝒉𝒕𝒆𝒅 𝒊𝒏𝒑𝒖𝒕𝒔 = 𝑺𝒖𝒎 𝒐𝒇 𝒘𝒆𝒊𝒈𝒉𝒕𝒆𝒅 𝒓𝒆𝒔𝒑𝒐𝒏𝒔𝒆

If a system does not satisfy this condition, then it is non-linear.

For discrete system, the condition for linearity is,

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

Time-Invariant System: A system is time invariant if the behavior and characteristics of

the system are fixed over time. A system is time invariant if a time shift in the input signal
results in an identical time shift in the output signal. For example, a time-invariant system
should produce 𝑦(𝑡 − 𝑡0 ) as the output when 𝑥(𝑡 − 𝑡0 ) is the input. Mathematically, it can
be written as,

𝐻[𝑥(𝑡 − 𝑡0 )] = 𝑦(𝑡 − 𝑡0 )

56
Code:

Let a system is defined as,

ℎ(𝑛) = 3𝑒𝑥𝑝−3𝑛 , 𝑓𝑜𝑟 𝑛 = 0,1,2,3,4

For this, the input signal is,

𝑥1 (𝑛) = [0 1 2 3 0 0 0]

𝑥2 (𝑛) = [0 0 1 1 1 0 0]

And the weights are,

𝑎=2

𝑏=3

We have to check whether the system is linear and time-invariant or not.

Then code for linearity test will require the previous convint function to be implemented
which is,

function[y,ty]=convint(x,tx,h,th)
%Inputs:
%x is the input signal vector
%tx is the times of the samples in x
%h is the impulse response vector
%th is times of the samples in h
%outputs:
%y is the output signal vector,
%length(y)=length(x)+length(h)-1
%ty is the time of the samples in y
dt=tx(2)-tx(1);
y=conv(x,h)*dt;
ty=(tx(1)+th(1))+[0:(length(y)-1)]*dt;

And the code for time-invariant check will require both convint and shift function to be
implemented. The convint function is provided above. Code for shift function is,

function y = shift(x,k)
if(k>0)
y = [zeros(1,abs(k)) x(1:end-abs(k))];
else
y = [x((abs(k)+1):end) zeros(1,abs(k))];
end
end

57
Then the code for LTI test is,

clc; clear;
x1 = [0 1 2 3 0 0 0];
x2 = [0 0 1 1 1 0 0];
a = 2;
b = 3;
len = length(x1);
tx=0:1:len-1;
th=0:1:4;
h =3* exp(-3*th);
[y1 ty1]=convint(x1,tx,h,th);
[y2 ty2]=convint(x2,tx,h,th);
LHS = a*y1+b*y2;

x3 = a*x1 + b*x2;
[y3 ty3]=convint(x3,tx,h,th);
RHS=y3;

LHS=round(LHS,4);
RHS=round(RHS,4);

subplot(2,2,1);
stem(tx,x1); % plot input-1
xlabel('Time');
ylabel('Amplitude');
title('First input sequence');

subplot(2,2,2);
stem(tx,x2); % plot input-2
xlabel('Time');
ylabel('Amplitude');
title('Second input sequence');
subplot(2,2,3);
stem(ty1,LHS); grid % plot LHS
xlabel('Time');
ylabel('Amplitude');
title(' LHS');

subplot(2,2,4);
stem(ty3,RHS); grid % plot RHS
xlabel('Time');
ylabel('Amplitude');
title(' RHS');

if(LHS == RHS)
disp('System is linear.');
else
disp('System is non-linear.');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%time-invariantcheck%%%%%%%%%%%%%%%%%%%
yshifted_output = shift(y1,2);

xshifted_input = shift(x1,2);

58
[yshifted_output2,ty4] = convint(xshifted_input,tx,h,th);

yshifted_output =round(yshifted_output,4);
yshifted_output2=round(yshifted_output2,4);

figure(2)
subplot(2,2,1);
stem(tx,x1); % plot input
title('Input sequence');

subplot(2,2,2);
stem(ty1,y1); % plot output
title('Output sequence');

subplot(2,2,3);
stem(ty1,yshifted_output); % plot only shifted output
title('Only Shifting the Output');
subplot(2,2,4);
stem(ty4,yshifted_output2); % plot output after shifting input
title('Output for Shifted version of input');

if(yshifted_output==yshifted_output2)
disp('System is Time Invariant.');
else
disp('System is not Time Invariant.');
end

59
Its output is,

Fig 4.2: Input and output signal for linearity test of this system.

Fig 4.3: Input and output signal for time invariance check test of this system.

60
Also,

Fig 4.4: Output of decision of the previous code.

Code:

Let a system is defined as,

𝑦(𝑛) = [𝑥(𝑛)]2

For this, the input signal is,

𝑥1 (𝑛) = [0 1 2 3 0 0 0]

𝑥2 (𝑛) = [0 0 1 1 1 0 0]

And the weights are,

𝑎=2

𝑏=3

We have to check whether the system is linear and time-invariant or not.

Then the code for it will be,

clc; clear;
x1 = [0 1 2 3 0 0 0];
x2 = [0 0 1 1 1 0 0];
a = 2;
b = 3;
len = length(x1);
tx=0:1:len-1;
y1= power(x1,2);
y2= power(x2,2);
LHS = a*y1+b*y2;

x3 = a*x1 + b*x2;
y3= power(x3,2);
RHS=y3;

LHS=round(LHS,4);
RHS=round(RHS,4);
ty=0:(length(y1)-1);
subplot(2,2,1);

61
stem(tx,x1); % plot input-1
xlabel('Time');
ylabel('Amplitude');
title('First input sequence');

subplot(2,2,2);
stem(tx,x2); % plot input-2
xlabel('Time');
ylabel('Amplitude');
title('Second input sequence');
subplot(2,2,3);
stem(ty,LHS); grid % plot LHS
xlabel('Time');
ylabel('Amplitude');
title(' LHS');

subplot(2,2,4);
stem(ty,RHS); grid % plot RHS
xlabel('Time');
ylabel('Amplitude');
title(' RHS');

if(LHS == RHS)
disp('System is linear.');
else
disp('System is non-linear.');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%time-
invariantcheck%%%%%%%%%%%%%%%%%%%%%%
yshifted_output = shift(y1,2);

xshifted_input = shift(x1,2);
yshifted_output2= power(xshifted_input,2);

yshifted_output=round(yshifted_output,4);
yshifted_output2=round(yshifted_output2,4);

figure(2)
subplot(2,2,1);
stem(tx,x1); % plot input
title('Input sequence');

subplot(2,2,2);
stem(ty,y1); % plot output
title('Output sequence');

subplot(2,2,3);
stem(ty,yshifted_output); % plot only shifted output
title('Only Shifting the Output');
subplot(2,2,4);
stem(ty,yshifted_output2); % plot output after shifting input
title('Output for Shifted version of input');

62
if(yshifted_output==yshifted_output2)
disp('System is Time Invariant.');
else
disp('System is not Time Invariant.');
end

Its output is,

Fig 4.5: Input and output signal for linearity test of this system.

63
 Fig 4.6: Input and output signal for time invariance check test of this system.

Fig 4.7: Output of decision of the previous code.

Verify whether the following systems are LTI or not,

 [𝑦(𝑛)]2 = 𝑛[𝑥(𝑛)]2  Non-linear and not Time Invariant.

 [𝑦(𝑛)]2 = 𝑛[𝑥(𝑛)] Linear and not Time Invariant.

64
Laplace Transform: The Laplace transform is a generalization of the Fourier transform
of a continuous time signal. The Laplace transform converges for signals such as unit ramp
signal, parabolic signal etc. for which the Fourier transform does not converge due to a
convergent factor 𝑒 −𝜎𝑡 . Hence, the Laplace transform is a useful tool in the analysis and
design of continuous time systems. The Laplace transform of a signal, 𝑓(𝑡)is,
∞
𝐹(𝑠) = ∫ 𝑓(𝑡)𝑒 −𝑠𝑡 𝑑𝑡
0

Where,

𝑠 = 𝜎 + 𝑗𝜔

Thus Laplace transform converts the time domain signal 𝑓(𝑡) into frequency domain
signal 𝐹(𝑠).

Existence of Laplace Transforms: For the existence of Laplace transform, the integral form
∞
𝐹(𝑠) = ∫0 𝑓(𝑡)𝑒 −𝑠𝑡 𝑑𝑡 must converge i.e. the Laplace transform is said to exist if the
magnitude of the transform is finite, that is, |𝐹(𝑠)| ≪ ∞. This limits the variable
𝑠 = 𝜎 + 𝑗𝜔 to a part of complex plane called the Region of Convergence (ROC). ROC is
required for inverse Laplace transform. For example, if

𝑓(𝑡) = 𝑒 3𝑡

1
𝑇ℎ𝑒𝑛, 𝐹(𝑠) =
𝑠−3

Here, the Laplace transform is defined only for,

𝑅𝑠 (𝑠) > 3

𝑖. 𝑒. 𝜎 > 3

65
 1
Fig 4.7: ROC for 𝐹(𝑠) = 𝑠−3

Inverse Laplace Transform: The inverse Laplace transform is defined as,

1 𝜎0 +𝑗𝜔
𝑓(𝑡) = ∫ 𝐹(𝑠)𝑒 𝑠𝑡 𝑑𝑠
2𝜋𝑗 𝜎0 −𝑗𝜔

MATLAB command for Laplace and inverse Laplace transform is,

𝑙𝑎𝑝𝑙𝑎𝑐𝑒()

𝑖𝑙𝑎𝑝𝑙𝑎𝑐𝑒()

Code: Code to find the Laplace transform of

𝑓(𝑡) = 𝑒 −𝑎𝑡

And to find the inverse Laplace transform of 𝐹(𝑠) is,

clc;
clear all;
syms t s a w % creates symbolic variables

f = exp(-a*t) % from Table 3.1

F = laplace(f) % computes the Laplace transform
simplify(F) % simplify the transformed equation
pretty(ans) % prints symbolic output as typeset mathematics.

f = ilaplace(ans)

66
Output is,

Fig 4.8: Output of above code for Laplace and inverse Laplace transform

Perform Laplace transform and thereby inverse Laplace transform for the below
functions:

 𝑓 = sin 𝜔𝑡
 𝑓 = 𝑡 sin 𝜔𝑡
 𝑓 = 𝑡 3 + 3𝑡 2 − 6𝑡 + 4
 𝑓 = 𝑐𝑜𝑠 3 𝜔𝑡

67
Z-Transform: The z-transform of an arbitrary discrete-time signal 𝑥(𝑛) is defined as:
𝑛=∞

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

𝑊ℎ𝑒𝑟𝑒, 𝑧 𝑖𝑠 𝑎 𝑐𝑜𝑚𝑝𝑙𝑒𝑥 𝑣𝑎𝑟𝑖𝑎𝑏𝑙𝑒.

It is two sided Z-transform. The expression for one-sided Z-transform is,

𝑛=∞

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

Inverse Z-transform: It is performed to derive the original time domain signal which is
expressed as,

𝑥(𝑛) = 𝑍 −1 [𝑋(𝑧)]

MATLAB command for Z and inverse Z-transform is,

𝑧𝑡𝑎𝑛𝑠()

𝑖𝑧𝑡𝑟𝑎𝑛𝑠()

Code: Code to find the Z-transform of

𝑓(𝑛) = 𝑒 −𝑎𝑛

And to find the inverse Z-transform of 𝐹(𝑧)

clc;
syms n z a w

f = exp(-a*n);
Z_tranz = ztrans(f) % computes the z transform
SZ_tranz = simplify(Z_tranz)
pretty(SZ_tranz)
f = iztrans(Z_tranz); % computes inverse z transform
pretty(ans)

68
Output is,

Fig 4.8: Output of above code for Z and inverse Z transform

Perform Z transform and thereby inverse Z transform for the below functions:

 𝑓(𝑛) = 𝑒 𝑎𝑛
 𝑓(𝑛) = 𝑛𝑒 −𝑎𝑛
 𝑓(𝑛) = sin 𝜔𝑛
 𝑓(𝑛) = cos 𝜔𝑛
 𝑓(𝑛) = 𝑒 −𝑎𝑛 sin 𝜔𝑛

Perform inverse Z transform for the below functions:

−4𝑧 2 +8𝑧
 𝐹 = 𝑧 2 +6𝑧+8
𝑧3
 𝐹 = (𝑧+1)(𝑧−1)2
𝑧
 𝐹 = 3𝑧 2 −4𝑧+1
2𝑧 3 +3𝑧 2
 𝐹 = (𝑧+1)∗(𝑧+0.5)∗(𝑧−2.5)

69
Stability of a system: If Poles lies within the unit circle of the zplane, then the system is
stable.

Code: Determine whether the following system is stable or not which has Z-transform for
a particular input as,

𝑧 2 − 2.5𝑧 + 1
𝑋(𝑧) =
(𝑧 2 − 0.75𝑧 + 0.075)

clc
num=[1 -2.5 1]; den=[1 -.75 0.075];
z = tf(num, den, 0.02) %tf ---> transfer function; Sample time Ts
= 0.02
[z, p, k]=tf2zpk(num, den) % finds zeros, poles, and gains of
transfer function.
[r, p, k]=residuez(num, den) % residuez converts a system,
expressed as the ratio of two polynomials, to partial fraction
expansion; r =residues, p = poles, k= direct terms.

zplane(num, den) % plots the zeros and poles in zplane

if abs(p)<1
disp('The System Is Stable.');
else
disp('The System Is Unstable.');
end

Its output is,

Fig 4.9: Poles-Zeros plot of the above system’s transfer function

70
And,

Fig 4.10: Decision if the system is stable or not

71
4.3. Experimental Task:

(1 Laplace transform; 2 inverse z-transform; 3 Stability test)

1. f = t*sin(a*t)
2. F = (z)/((z-1)*(z-2))
3. F = (2z^3 + 3z^2 - 3.2z - 1.6)/(2z^3 + 2.5z^2 + 0.25z - 0.25);

1. f = (1-exp(t))/t
2. F = z^2/((z-1)*(z-0.2))
3. F = (2z^3 - 3.5z^2 - 2.5z + 3)/(2z^3 - 2z^2 - 0.5z + 0.5);

72
Experiment 05: Study on Continuous and Discrete State Space Model.
5.1. Objectives:
1. To implement continuous and discrete state space model and observe its step response,
magnitude and phase response.
2. To check the stability, controllability and observability of the system.

5.2. Theory:

A state-space model is a first-order differential equation-based mathematical

representation of a physical system that includes input, output, and state variables. The
values of the output variables are determined by the state variables.

For time invariant system,

̇ + 𝐵𝑢(𝑡)
𝑥(𝑡) = 𝐴𝑥(𝑡)

𝑦(𝑡) = 𝐶𝑥(𝑡) + 𝐷𝑢(𝑡)

𝑊ℎ𝑒𝑟𝑒,

𝐴 = (𝑁𝑥 ∗ 𝑁𝑥 ) 𝑚𝑎𝑡𝑟𝑖𝑥 𝑘𝑛𝑜𝑤𝑛 𝑎𝑠 𝑒𝑣𝑜𝑙𝑢𝑡𝑖𝑜𝑛 𝑚𝑎𝑡𝑟𝑖𝑥

𝐵 = (𝑁𝑥 ∗ 𝑁𝑢 ) 𝑚𝑎𝑡𝑟𝑖𝑥 𝑘𝑛𝑜𝑤𝑛 𝑎𝑠 𝑐𝑜𝑛𝑡𝑟𝑜𝑙 𝑚𝑎𝑡𝑟𝑖𝑥

𝐶 = (𝑁𝑦 ∗ 𝑁𝑥 ) 𝑚𝑎𝑡𝑟𝑖𝑥 𝑘𝑛𝑜𝑤𝑛 𝑎𝑠 𝑜𝑏𝑠𝑒𝑟𝑣𝑎𝑡𝑖𝑜𝑛 𝑚𝑎𝑡𝑟𝑖𝑥

𝐴 = (𝑁𝑦 ∗ 𝑁𝑢 ) 𝑚𝑎𝑡𝑟𝑖𝑥 𝑘𝑛𝑜𝑤𝑛 𝑎𝑠 𝑑𝑖𝑟𝑒𝑐𝑡 𝑡𝑟𝑎𝑛𝑠𝑚𝑖𝑠𝑠𝑖𝑜𝑛 𝑚𝑎𝑡𝑟𝑖𝑥

For implementing in MATLAB,

 𝑠𝑠(𝐴, 𝐵, 𝐶, 𝐷) is used for implementation of continuous state space model.

 𝑠𝑠(𝐴, 𝐵, 𝐶, 𝐷, 𝑡𝑠) is used for implementation of discrete state space model. To leave the
sample time unspecified, 𝑡𝑠 is set to -1.

73
Code for Continuous State Space Model:

Implement the state space model of the following system,

Fig 5.1: Circuit Diagram of a RLC circuit

Where it is given that,

1
0 −
𝐴=[ 𝐶]
1 𝑅
−
𝐿 𝐿
1/𝐶
𝐵=[ ]
0
𝐶 = [0 𝑅]
𝐷=0
Let’s assume, 𝑅 = 100Ω, 𝐶 = 1𝜇𝐹, 𝐿 = 0.22𝜇𝐻

Then the code will be,

clc;
clear;

cp=1*10^(-9);
L=0.22*10^(-9);
res=100;

A=[0 -1/cp;
1/L -res/L];
B=[1/cp; 0];
C=[0 res];
D=0;
sys=ss(A,B,C,D)

figure(1)
step(sys)
figure(2)

74
bode(sys)

stb=isstable(sys);
if stb==1
disp('Stable')
else
disp('Not Stable')
end

S = ctrb(sys);
Co=det(S);
if Co==0
disp('Not Controllable')
else
disp('Controllable')
end

V = obsv(sys);
ob=det(V);
if ob==0
disp('Not Observable')
else
disp('Observable')
end

Its output is,

75
 Fig 5.2: Step Response of this system

Fig 5.3: Magnitude and Phase response of the system

76
Fig 5.4: Stability, Controllability and Observability check of this system

77
Code for Discrete State Space Model:

Implement the discrete state space model when the following state matrices are given:

0 1 0 0
0 0 1 0
𝐴=[ ]
0 0 0 1
−3 −6 −5 −4

0
0
𝐵=[ ]
0
1

𝐶 = [1 0 0 0]

𝐷=0

The code for this will be,

clc;
clear;

A=[0 1 0 0;
0 0 1 0;
0 0 0 1;
-3 -6 -5 -4];
B=[0; 0; 0 ; 1];
C=[1 0 0 0];
D=0;
sys=ss(A,B,C,D,-1)

figure(1)
step(sys)
figure(2)
bode(sys)

stb=isstable(sys);
if stb==1
disp('Stable')
else
disp('Not Stable')
end

S = ctrb(sys);
Co=det(S);
if Co==0

78
 disp('Not Controllable')
else
disp('Controllable')
end

V = obsv(sys);
ob=det(V);
if ob==0
disp('Not Observable')
else
disp('Observable')
end

Its output is,

Fig 5.5: Step Response of this system

79
 Fig 5.6: Magnitude and Phase response of the system

Fig 5.7: Stability, Controllability and Observability check of this system

80