0% found this document useful (0 votes)
53 views

DSP Lab Manual Final Final

This document is the lab manual for the Digital Signal Processing course at NED University of Engineering & Technology. It contains instructions for 12 experiments related to key concepts in digital signal processing, including sampling, quantization, discrete-time signals and systems, the discrete Fourier transform, and digital filter design. The experiments are designed to help students verify analytical techniques and achieve the program learning outcomes for discrete-time systems and signals.

Uploaded by

Taskeen Junaid
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as DOC, PDF, TXT or read online on Scribd
0% found this document useful (0 votes)
53 views

DSP Lab Manual Final Final

This document is the lab manual for the Digital Signal Processing course at NED University of Engineering & Technology. It contains instructions for 12 experiments related to key concepts in digital signal processing, including sampling, quantization, discrete-time signals and systems, the discrete Fourier transform, and digital filter design. The experiments are designed to help students verify analytical techniques and achieve the program learning outcomes for discrete-time systems and signals.

Uploaded by

Taskeen Junaid
Copyright
© © All Rights Reserved
We take content rights seriously. If you suspect this is your content, claim it here.
Available Formats
Download as DOC, PDF, TXT or read online on Scribd
You are on page 1/ 89

NED University of Engineering &

Technology Department of Electrical


Engineering

LAB MANUAL
For the course

Digital Signal Processing


(EE-394) For T.E.(EE)

Instructor name:__________________________________
Student name:____________________________________
Roll no: Batch:___________________
Semester: Year:________________
LAB MANUAL
For the course

Digital Signal Processing


(EE-394) For T.E.(EE)

Content Revision Team:


Mr. Muhammad Omar & Mr. Nabeel Fayyaz
Last Revision Date: 29-12-2020

Approved By

The Board of Studies of Department of Electrical Engineering

____________________ ____________________

____________________ ___________________

____________________ ____________________
CONTENTS
Psychomotor Level: 4
CLO: Experimentally verify the analytical and design techniques developed for discrete time systems and signals.
[PLO (5)]
S.No. Date Title of Experiment Total Marks Signature

1 To study the effects of Sampling in Discrete Time


Signals.

2 To study the effects of Quantization in Discrete Time


Discrete Valued Signals.

3 To study and verify Discrete-Time convolution and its


properties.

4 To study Discrete-Time correlation.

5 The Discrete Fourier Transform as a Linear


Transformation

6 Studying Discrete Fourier Transform using an audio


signal example

7 DFT or FFT leakage and its solution.

8 s-plane and plot impulse and frequency response for


different pole zero location in s-plane.
z-plane and plot impulse and frequency response
9 for different pole zero location in z-plane.

10 Designing FIR filters with windowing and designing


IIR filters using FDA tool

11 Open Ended Lab1 (due after lab 2)

12 Open Ended Lab2(due after lab 6)


LAB SESSION 01
OBJECTIVE:
To study the relationship between discrete-time and continuous time signals by examining
sampling and aliasing.
THEORY:
Signals are physical quantities that carry information in their patterns of variation. Continuous-
time signals are continuous functions of time, while discrete-time signals are sequences of
numbers. If the values of a sequence are chosen from a finite set of numbers, the sequence is
known as a digital signal. Continuous-time, continuous-amplitude signals are also known as
analog signals.
Analog phenomenon is continuous – like a human speech of speaker, or a continuously rotating
disc attached to the shaft of motor etc. With analog phenomena, there is no clear separation
between one point and the next; in fact, between any two points, an infinite number of other
points exist. Discrete phenomenon, on the other hand, is clearly separated. There's a point (in
time or space), and then there's a neighboring point, and there's nothing between the two.
Signal processing is concerned with the acquisition, representation, manipulation, transformation,
and extraction of information from signals. In analog signal processing these operations are
implemented using analog electronic circuits. Converting the continuous phenomena of images,
sound, and motion into a discrete representation that can be handled by a computer is called analog-
to-digital conversion. Digital signal processing involves the conversion of analog signals into digital,
processing the obtained sequence of finite precision numbers using a digital signal processor or
general purpose computer, and, if necessary, converting the resulting sequence back into analog
form. When stored in a digital computer, the numbers are held in memory locations, so they would
be indexed by memory address. Regardless of the medium (either sound or an image), analog-to-
digital conversion requires the same two steps: Sampling and Quantization.
Sampling: This operation chooses discrete (finite) points at which to measure a continuous
phenomenon (which we will also call a signal). In the case of sound, the sample points are
evenly separated in time. In the case of images, the sample points are evenly separated in space.
Sampling Rate: The number of samples taken per unit time or unit space is called the sampling
rate. The frequency of sampled/discrete phenomenon (signal) can be calculated as
fd = F /Fs (cycles/sec )/(samples/sec) = cycles/ samples
Where, F = Frequency of analog or continuous phenomenon (signal). [Unit: cycles/sec]
Fs = Sampling frequency or sampling rate [Unit: samples/sec]
fd = Frequency of Discrete phenomenon (signal). [Unit: cycles/sample]
Sampling Theorem: A continuous time phenomenon or signal like x(t) can be reconstructed
exactly from its samples x(n) = x(nTs), if the samples are taken at a rate Fs = 1/Ts
that is greater than twice the frequency of the signal being sampled i.e. Fs ≥ 2 ∗ F. Mathematically,
Aliasing: A common problem that arises when sampling a continuous signal is aliasing, where a
sampled signal has replications of its sinusoidal components which can interfere with other
components. It is an effect that causes two discrete time signals to become indistinct due to
improper sampling (fd>1/2 cycles/sample).
PROCEDURE:
1. Simulate and plot two CT signals of 10 Hz and 110 Hz for 0 < t < 0.2 secs.
2. Sample at Fs = 100 Hz and plot them in discrete form.
3. Observe and note the aliasing effects.
4. Explore and learn.
STEPS:
1. Make a folder at desktop and name it as your current directory within MATLAB.
2. Open M-file editor and type the following code:
clear all;
close all;
clc;
F1 = 10;
F2 = 110;
Fs = 100;
Ts = 1/Fs;
t = [0 : 0.0005 : 0.2];
x1t = cos(2*pi*F1*t);
x2t = cos(2*pi*F2*t);
figure,
plot(t,x1t,t,x2t, 'LineWidth',2);
xlabel('continuous time (sec)');
ylabel('amplitude');
xlim([0 0.1]);
grid on;
legend('10Hz','110Hz');
title('Two CTCV(analog) sinusoids');
3. Save the file as P011.m in your current directory and ‘run’ it, either using F5 key or
writing the file name at the command window.
(Check for the correctness of the time periods of both sinusoids.)
Now add the following bit of code at the bottom of your P011.m file and save.

nTs = [0 :Ts : 0.2]; ylabel('amp');


n = [1 : length(nTs)-1 ]; xlim([0 0.1]);
x1n = cos(2*pi*F1*nTs); subplot(2,1,2)
x2n = cos(2*pi*F2*nTs); stem(nTs,x2n,'LineWidth',2);
figure, grid on;
subplot(2,1,1), title('110Hz sampled')
stem(nTs,x1n,'LineWidth',2); xlabel('discrete time(sec)');
grid on; ylabel('amp');
xlabel('discrete time (sec)'); xlim([0 0.1]);

1. Before hitting the ‘run’, just try to understand what the code is doing and try to link it
with what we have studied in classes regarding concepts of frequency for DT signals.

2. Now ‘run’ the file and observe both plots.


To see what is really happening, type the following code at the bottom of your existing
P011.m file and run again.
figure, xlabel('time (sec)');
plot(t,x1t,t,x2t); ylabel('amplitude');
hold; xlim([0 0.05]);
stem(nTs,x1n,'r','LineWidth',2); legend('10Hz','110Hz');
3. Observe the plots.
RESULT:
Explain (write) in your own words the cause and effects of what you just saw.

LAB TASKS:
1. Consider the following CT signal:
x(t) = sin (2 pi F0 t). The sampled version will be:
x(n) = sin (2 pi F0/Fs n),
where n is a set of integers and sampling interval Ts=1/Fs.
Plot the signal x(n) for n = 0 to 99
for Fs = 5 kHz and F1 = 0.5, 2, 3 and 4.5 kHz. Explain the similarities and differences among
various plots.
clear all;
close all;
clc;
F1=500; %2000; %3000; %4500; %( here we are changing frequency of
analog signal keeping sampling frequency same to observe their plots
on same sampling frequency)
Fs=5000;
Ts=1/Fs
nTs=[0: Ts :99]
x1n=sin(2*pi*F1*nTs)
figure
stem(nTs,x1n,'filled','LineWidth',1)
grid on
title('500 Hz DTCV (sampled) signal')
xlabel('discrete time (sec)')
ylabel('amplitude')
xlim([0 0.01])
2. Generate a tone in MATLAB with varying frequency f = 1000,2000,3000,4000, 5000, 6000, 8000,
9000, 25000,-1000,-2000,-3000 Hz with Fs = 8000 samples/sec.
Listen to the tones, and observe at Sounds like what
frequency? Also Specify whether Aliasing is happening or not.
CODE FOR TASK 2:
Fs=8000;
f1=1000; f2=2000; f3=3000; f4=4000; f5=5000; f6=6000; f7=8000; f8=9000; f9=25000; f10=-
1000; f12=-2000; f13=-3000;
t= [0:1/Fs:5];
x= sin(2*pi*f2*t);
sound(x)
OBSERVATION:
Yes aliasing is happening, after listening the tones at different analog frequencies we can
say that.
f3 is an aliasing effect of f5 and f13. f1 is an aliasing effect of f8,f9 and f10.
f2 is an aliasing effect of f6 and f12.

3. Record a sentence in your voice.(you may use Simulink /audacity to record).Change Fs =44100,
22050, 11025, 8192, 4096 , 2048 , 1024
and observe
a) Voice quality during playback [Excellent/Good/OK/Bad]
b) File size in kilobytes
c) Aliasing happening or not?

OBSERVATION:
Voice quality is turning bad as we decrease the sampling frequency Yes aliasing is happening, after
decreasing sampling one by one frequency we can say that signal is losing its information.

OUTCOME BASED EDUCATION


LAB RUBRICS
Digital Signal Processing (EE-394)
Student Name: Roll No: Section:

Excellent Good Average Poor Comments


100% 75% 50% 0%
The background The concepts The theoretical There is no
knowledge of
Clarity of theory is fully are not knowledge or
the background
Concepts grasped. completely concepts are
theory/concepts.
cleared. inappropriate.
The logic of
MATLAB the
program/Simulink program/model
The The
model is properly is correct but
Obtaining program/model program/model
simulated to there is some
plot/results has some is not correct.
obtain the desired syntax error to
major errors.
response. obtain the
results

Lab task Late


Submission submission of
submitted in time.
Lab Task
LAB SESSION 02
OBJECTIVE:
To observe the quantization effects on sampled signals and to understand how quantization leads to
quantization error. In this lab, we will investigate the influence of the number of quantization levels
on the quality of digitized signal. Method of selection of ADC is also a part of this lab session.
THEORY:
Everything stored on a computer is discrete time discrete valued signal. Because computer has finite number of
registers and each register is a finite length register. We take too many samples to give the ‘effect’ of continuous
time signals. But actually they are discrete time. We also take very fine resolution of amplitude axis to give the
effect of continuous valued signal but due to finite word length of the computer register, the stored variables are
already quantized. This lab aims to explain the quantization effects in a computer.
Regardless of the medium (audio or image), the digitization of real world analog signal usually involves
two stages: sampling, i.e. the measurement of signal at discretely spaced time intervals, and quantization,
i.e. the transformation of the measurements (amplitudes) into finite -precision numbers (allowed discrete
levels), such that they can be represented in computer memory. Quantization is a matter of representing
the amplitude of individual samples as integers expressed in binary. The fact that integers are used forces
the samples to be measured in a finite number of bits ( discrete levels). The range of the integers possible is
determined by the bit depth, the number of bits used per sample. The bit depth limits the precision with
which each sample can be represented.
Bit Depth:
Within digital hardware, numbers are represented by binary digits known as bits—in fact, the term bit
originated from the words Binary digit. A single bit can be in only one of two possible states: either a one
or a zero. When samples are taken, the amplitude at that moment in time must be converted to integers in
binary representation. The number of bits used to represent each sample, called the bit depth
(bits/sample) or sample size, determines the precision with which the sample amplitudes can be
represented. Each bit in a binary number holds either a 1 or a 0. In digital sound, bit depth affects how
much you have to round off the amplitude of the wave when it is sampled at various points in time.
The number of different values that can be represented with b-bit is 2 b .The largest decimal number that
can be represented with an b-bit binary number is 2 b - 1. For example, the decimal values that can be
represented with an 8-bit binary number range from 0 to 255, so there are 256 different values (levels of
ADC). A bit depth of 8 allows 28=256 different discrete levels at which samples can be approximated or
recorded. Eight bits together constitute one byte. A bit depth of 16 allows 2 16 = 65,536 discrete levels,
which in turn provides much higher precision than a bit depth of 8.
The number of bits in a data word is a key consideration. The more bits used in the word, the better the
resolution of the number, and the larger the maximum value that can be represented. Some computers use
64-bit words. Now, 264 is approximately equal to 1.8 x 1019—that's a pretty large number. So large, in fact,
that if we started incrementing a 64-bit counter once per second at the beginning of the universe (≈20
billion years ago), the most significant four bits of this counter would still be all zeros today.
To simplify the explanation, take an example of ADC with a bit depth of 3, 2 3 = 8 quantization levels
ranging from -4 to 3 are possible in signed magnitude representation. For bipolar ADCs (or signed
magnitude representation), by convention, half of the quantization levels are below the horizontal axis (that
is 21, of the quantization levels). One level is the horizontal axis itself (level 0), and 2b-1 − 1levels are
above the horizontal axis. Note that since one bit is used for the signed bit (in 2-complementformat), the
largest magnitude corresponds to 2^(b -1 ). (not 2b). When a sound is sampled, each sample must be
scaled to one of the 8 discrete levels. However, the samples in reality might not fall neatly onto these
levels. They have to be rounded up or down by some consistent convention.
QUANTIZATION ERROR:
The samples, which are taken at evenly-spaced points in time, can take on the values only at the discrete
quantization levels to store on our computer. Therefore quantization leads to a loss in the signal quality,
because it introduces a “Quantization error”. Quantization error is sometimes referred to as
'"Quantization noise". Noise can be broadly defined as part of an audio signal that isn’t supposed to be
there. However, some sources would argue that a better term for quantization error is " distortion",
defining distortion as an unwanted part of an audio signal that is related to the true signal.
The difference between the quantized samples and the original samples constitutes quantization error or
rounding error (if round-off method is used). Xe(n) = Xq(n) − x(n). The lower the bit depth, the more
values potentially must be approximated (rounded), resulting in greater quantization error.
To calculate the required bit depth of ADC i.e. bits/sample, there are two important points which we must
have to consider:
a) How much noise is already present in the analog signal?
b) How much more noise can be tolerated in the digital signal?
Signal-to -noise-ratio- SNR (of analog signal)
Before looking at SNR specifically in the context of digital imaging and sound, let's co nsider the general
definition. Signal-to-noise ratio can generally be defined as the ratio of the meaningful content of a signal
versus the associated background noise.
SNR = 10log10 (Px /Pe )
Where, Px and Pe are average power of the analog signal and noise signal respectively.
A signal is any type of communication – something a person says to you, a digital signal sending an
image file across a network, a message posted on an electronic bulletin board, a piece of audio being
played from a cassette, etc. The noise is the part of the message that is not meaningful; in fact, it gets in
the way of the message intended in the communication. You could use the term signal-to-noise ratio to
describe communications in your everyday life. If you know someone who talks a lot but doesn't really
convey a lot of meaning, you could say that he or she has a low signal - to-noise ratio. Web-based bulletin
board and chat groups are sometimes described as having a low SNR – there may be quite a few postings,
but very much meaningful content. In these first two examples, the noise consists of all the empty "filler"
words. In the case of a digital signal sent across a network, the noise is the electronic degradation of the
signal. On a piece of audio played from cassette, the noise could be caused by damage to the tape or
mechanical imperfections in the cassette player.

In analog data communication (analog signals), the signal -to-noise ratio is defined as the ratio of the
average power in the signal versus the power in the noise level. In this context, think of a signal being
sent over a network connection compared to the extent to which the signal is corrupted. This is related to
the general usage of the term described above. This usage of the term SNR applies to analog signals.

SIGNAL-TO-QUANTIZATION-NOISE-RATIO-SQNR (OFADC):
Using finite word lengths prevents us from representing values with infinite precision, increases the
background noise in our spectral estimation techniques etc. The amount of error implicit in a chosen bit
depth can be measured in terms of the signal-to-noise ratio (SNR).
For a digitized image or sound, the signal-to-noise ratio is defined as the ratio of the maximum
sample value versus the maximum quantization error. In this usage of the term, the ratio depends on the bit
depth chosen for the signal. Any signal encoded with a given bit depth will have the same ratio. This can
also be called signal-toquantization-noise ratio (SQNR), but you should be aware that in many sources the
term signal-to-noise ratio is used with this meaning as well. (Henceforth, we'll use the term SQNR to
distinguish this measurement from SNR.)
Practical A/D converters are constrained to have binary output words of finite length. Commercial A/D
converters are categorized by their output word lengths, which are normally in the range from 8 to 16 bits.
There is no infinite bit ADC. So whenever we will digitize our signal, we will always have a quantization
error. Quantization error represents the quality of quantization process but the total error may also turn
out to be zero, so signal-toquantization-noise-ratio (SQNR) is used to describe the quality of quantization
process and it can be defined as
SQNRA/D = 10 log 10(Px /Pe)
Where, Px~and Pe are average power of the DTCV (sampled) signal and quantization error signal
respectively.

PROCEDURE:
1. Simulate a DTCV sinusoid of 1/50 cycles/sample with length of the signal be 500.
2. Choose the no. of significant digits for round-off and apply to the signal generated above.
3. Compute the error signals and SQNR
4. Explore and observe.
STEPS:
1. Make a folder at desktop and name it as your current directory within MATLAB.
2. Open M-file editor and write the following code:
3. Save the file as P021.m in your current directory and run it.

clear all; SQNR = 10*log10(Px1/Pe1);


close all; disp(['The Signal to
clc; Quantization Noise Ratio is: '
fd1 = 1/50; num2str(SQNR) ' dB.' ]);
n = [0 : 499 ]; q=input('No. figure,
of Digits after decimal subplot(2,1,1);
points to be retained (0-9): plot(n,x1,n,x1q);
'); xlabel('indices')
x1 = cos(2*pi*fd1*n);
Px1 =
sum(abs(x1).^2)/length(x1);
x1q = round(x1*10^q)/10^q;
x1e = x1 -x1q;
Pe1 =
sum(abs(x1e).^2)/length(x1e);
ylabel('Amp'); xlim([0 49]);
ylim([-1.1 1.1]);
legend('DTCV','DTDV');
subplot(2,1,2);
plot(n,x1e);
xlabel('indices');
ylabel('Error');
xlim([0 49]);

No. of Digits after decimal points to


be retained(0-9):6

The signal to Quantization Noise Ratio is: 127.7491 dB

Explore and take notes.


Now modify the above code as follows and save as another file P022.m.
clear all;
close all;
clc;
fd1 = 1/50;
n = [0 : 499 ];
q = [0 : 10];
% No. of Digits after decimal
points to be retained for num =
1: for num = 1 : length(q)
x1 = cos(2*pi*fd1*n);
Px1 = sum(abs(x1).^2)/length(x1);
x1q=round(x1*10^q(num))/10^q(num);
x1e = x1q-x1;
Pe1= sum(abs(x1e).^2)/length(x1e);
SQNR(num) = 10*log10(Px1/Pe1); end

figure,
plot(q,SQNR);
xlabel('Significant
Digits'); ylabel('SQNR
(dB)'); xlim([q(1) q(end)]);
1. Before hitting the ‘run’, just try to understand what the code is doing and try to link it with
the previous code.
2. Now ‘run’ the file and observe the results.
RESULT:
Explain (write) in your own words the cause and effects of what you just saw.

LAB TASKS:
1. Effects of Quantization with variable precision levels
Simulate a DTCV sampled composite signal of 1=1/25 samples/sec and
2=1/50 samples/sec with length of the signal be 250 samples.
Take the desired number of significant digits from user as an input. Then choose the method
of Quantization (round-off, floor & ceil) and apply to the signal generated above. Compute the
quantization error signals and SQNR.
CODE FOR TASK:
fd1=1/25; Pe2r=sum(abs(x2r-x2).^2)/length(x2r-x2);
fd2=1/50; Pe2c=sum(abs(x2c-x2).^2)/length(x2c-x2);
n=[0:249]; Pe2f=sum(abs(x2f-x2).^2)/length(x2f-x2);
x1=sin(2*pi*fd1*n); SQNRr1=10*log(Px1/Pe1r);
x2=sin(2*pi*fd2*n); SQNRc1=10*log(Px1/Pe1c);
Px1=(sum(abs(x1).^2))/length(x1); SQNRf1=10*log(Px1/Pe1f);
Px2=(sum(abs(x2).^2))/length(x2); SQNRr2=10*log(Px2/Pe2r);
N=input('Enter the desired number of SQNRc2=10*log(Px2/Pe2c);
significant figures:'); SQNRf2=10*log(Px2/Pe2f);
x1r=round(x1*10^N)/10^N ; disp(['signal to quantization round off noise
x1c=ceil(x1*10^N)/10^N ; ratio for x1 is: ' num2str(SQNRr1)]);
x1f=floor(x1*10^N)/10^N ; disp(['signal to quantization ceiling noise
x2r=round(x2*10^N)/10^N ; ratio for x1 is: ' num2str(SQNRc1)]);
x2c=ceil(x2*10^N)/10^N ; disp(['signal to quantization floor noise ratio
x2f=floor(x2*10^N)/10^N ; for x1 is: ' num2str(SQNRf1)]);
Pe1r=sum(abs(x1r-x1).^2)/length(x1r-x1); disp(['signal to quantization round off noise
Pe1c=sum(abs(x1c-x1).^2)/length(x1c-x1); ratio for x2 is: ' num2str(SQNRr2)]);
Pe1f=sum(abs(x1f-x1).^2)/length(x1f-x1); disp(['signal to quantization ceiling noise
ratio for x2 is: ' num2str(SQNRc2)]); subplot(2,3,6)
disp(['signal to quantization floor noise ratio stem(n,x1f-x1);
for x2 is: ' num2str(SQNRf2)]); xlabel('Error')
figure, ylabel('Amp')
suptitle('For fd=1/25') xlim([0 10])
subplot(2,3,1) ylim([-0.2 0.2])
stem(n,x1); grid on
hold; figure,
stem(n,x1r); suptitle('fd=1/50')
xlabel('Indices') subplot(2,3,1)
ylabel('Amp') stem(n,x2);
xlim([0 10]) hold;
ylim([-1 1]) stem(n,x2r);
grid on xlabel('Indices')
subplot(2,3,4) ylabel('Amp')
stem(n,x1r-x1); xlim([0 10])
xlabel('Error') ylim([-1 1])
ylabel('Amp') grid on
xlim([0 10]) subplot(2,3,4)
ylim([-0.2 0.2]) stem(n,x2r-x2);
grid on xlabel('Error')
subplot(2,3,2) ylabel('Amp')
stem(n,x1); xlim([0 10])
hold; ylim([-0.2 0.2])
stem(n,x1c); grid on
xlabel('Indices') subplot(2,3,2)
ylabel('Amp') stem(n,x2);
xlim([0 10]) hold;
ylim([-1 1]) stem(n,x2c);
grid on xlabel('Indices')
subplot(2,3,5) ylabel('Amp')
stem(n,x1c-x1); xlim([0 10])
xlabel('Error') ylim([-1 1])
ylabel('Amp') grid on
xlim([0 10]) subplot(2,3,5)
ylim([-0.2 0.2]) stem(n,x2c-x2);
grid on xlabel('Error')
subplot(2,3,3) ylabel('Amp')
stem(n,x1); xlim([0 10])
hold; ylim([-0.2 0.2])
stem(n,x1f); grid on
xlabel('Indices') subplot(2,3,3)
ylabel('Amp') stem(n,x2);
xlim([0 10]) hold;
ylim([-1 1]) stem(n,x2f);
grid on xlabel('Indices')
ylabel('Amp')
xlim([0 10]) xlabel('Error')
ylim([-1 1]) ylabel('Amp')
grid on xlim([0 10])
subplot(2,3,6) ylim([-0.2 0.2])
stem(n,x2f-x2); grid on

Enter the desired number of significant figures:4


signal to quantization round off noise ratio for x1 is: 201.7196
signal to quantization ceiling noise ratio for x1 is: 188.2725
signal to quantization floor noise ratio for x1 is: 188.8919
signal to quantization round off noise ratio for x2 is: 201.7196
signal to quantization ceiling noise ratio for x2 is: 188.6395
signal to quantization floor noise ratio for x2 is: 188.515

2.Simple sinusoid quantized to various bits per sample:


Generate a 100 Hz sinusoid sampled at 10000 samples/sec and quantized at 1_bit/sample. Now
increase the bit depth for various numbers of bits per sample (2, 3, 4, 5, 6, 7, 8) and attach
plots. You can use two column format for plotting (but the diagrams should be visible).
Code for task:
clear all; x= sin(2*pi*fd*n);
close all; xq=round(x*10^q)/10^q;
clc; plot(n,xq);
fd= 1/100; xlabel('indices');
n= [0:499]; ylabel('Amp');
q=input('no. of bits per sample:') title('At N bit/sample');
3. Audio signal quantization to various bits per sample
Use your recorded voice in last session and quantize it at 1 bit /sample.
Change bit depth to 2,3,4 and then listen and take notes of your observations.
Decide no. of bits for audio until quality stops improving.
Observation:
As we increase the number of bits the sound becomes more clear and we are getting more data because
the more bits we use, the higher the resolution becomes.

OUTCOME BASED EDUCATION


LAB RUBRICS
Digital Signal Processing (EE-394)
Student Name: Roll No: Section:

Excellent Good Average Poor Comments


100% 75% 50% 0%
The background The concepts The theoretical There is no
knowledge of
Clarity of theory is fully are not knowledge or
the background
Concepts grasped. completely concepts are
theory/concepts.
cleared. inappropriate.
The logic of
MATLAB the
program/Simulink program/model
The The
model is properly is correct but
Obtaining program/model program/model
simulated to there is some
plot/results has some is not correct.
obtain the desired syntax error to
major errors.
response. obtain the
results

Lab task Late


Submission submission of
submitted in time.
Lab Task
LAB SESSION 03
OBJECTIVE:
To study impulse response, observe convolution technique in signal processing, and verify different properties
like causality, commutative, distributive and associative properties.

THEORY:

1. Convolution is given as :y(n) = x(n)*h(n) =

i.e.one can compute the output y(n) to a certain input x(n) when impulse response h(n) of that system is
known. Convolution holds commutative property.

2. The length of the resulting convolution sequence is N+M-1,where N and M are the lengths of two
convolved signals respectively.

3. In causal system, the outputs only depend on the past and/or present values of inputs and NOT on future
values. This means that the impulse response h(n) of a causal system will always exist only for n≥ 0.

PROCEDURE:

1. We have the impulse response of a system as h(n) = { 3,2,1,-2,1,0,-4,0,3}



2. For x(n)={1,-2,3,-4,3,2,1}

STEPS:
1. Make a folder at desktop and name it as your current directory within MATLAB.
2. Open M-file editor and write the following code:
clear all;
close all;
clc;
h = [3 2 1 -2 1 0 -4 0 3]; % impulse response
org_h = 1; % Sample number where origin exists
nh = [0 : length(h)-1]- org_h + 1;
x = [1 -2 3 -4 3 2 1]; % input sequence
org_x = 1; % Sample number where origin exists
nx = [0 : length(x)-1]- org_x + 1; y =
conv(h,x);

ny = [nh(1)+ nx(1) : nh(end)+nx(end)];


figure,

subplot(3,1,1),
stem(nh,h);
xlabel('Time index n');
ylabel('Amplitude');
xlim([nh(1)-1 nh(end)+1]);
title('Impulse Response h(n)');
grid;

subplot(3,1,2),
stem(nx,x);
xlabel('Time index n');
ylabel('Amplitude');
xlim([nx(1)-1 nx(end)+1]);
title('Input Signal x(n)');
grid;

subplot(3,1,3)
stem(ny,y);
xlabel('Time index n');
ylabel('Amplitude');
xlim([ny(1)-1 ny(end)+1]); title('Output
Obtained by Convolution'); grid;

1. Save the file as P031.m in your current directory and ‘run’ it.
2. Calculate the length of input signal (N) and impulse response (M) used in above task?
3. Calculate the length of the output sequence and verify the result with N+M-1
4. Try to learn, explore the code and make notes.
5. Now modify the above code such that h(n)= {3,2, 1, -2,1,0,-4,0,3}(origin is shifted) and check for
causality.
RESULT:

EXERCISE:

1. What will happen if we input x(n)={0,0,1,0,0} into the above system.



2. Can you prove the commutative property of the convolution?

3. Modify the code to prove Associative and Distributed properties of the convolution.

4. Convolve your recorded sound with drumloop.wav. Note your observation

a) Plot the output.


b) Listen the output
CODE AND OUPUT FOR TASK 1:
h = [3 2 1 -2 1 0 -4 0 3]; grid on
org_h=1; subplot(3,1,2)
nh=[0:length(h)-1]-org_h+1; stem(nx,x);
x=[0,0,1,0,0]; xlabel('Time index n');
org_x=3; ylabel('Amplitude');
nx=[0:length(x)-1]-org_x+1; title('Input response x(n)');
y=conv(h,x); xlim([nx(1)-1 nx(end)+1])
ny=[nh(1)+nx(1):nh(end)+nx(end)]; grid on
figure, subplot(3,1,3)
subplot(3,1,1) stem(ny,y);
stem(nh,h); xlabel('Time index n');
xlabel('Time index n'); ylabel('Amplitude');
ylabel('Amplitude'); title('Output of convolution y(n)');
title('Impulse response h(n)'); xlim([ny(1)-1 ny(end)+1])
xlim([nh(1)-1 nh(end)+1]) grid on
CODE AND OUTPUT FOR TASK 2:

h = [3 2 1 -2 1 0 -4 0 3]; xlabel('Time index n');


org_h=1 ylabel('Amplitude');
nh=[0:length(h)-1]-org_h+1; title('Output of convolution y(n)=h(n)*x(n)');
x=[1 -2 3 -4 3 2 1]; xlim([ny(1)-1 ny(end)+1])
org_x=1; grid on
nx=[0:length(x)-1]-org_x+1; subplot(2,1,2)
y=conv(h,x); stem(ny,y);
ny=[nh(1)+nx(1):nh(end)+nx(end)]; xlabel('Time index n');
y1=conv(x,h) ylabel('Amplitude');
ny1=[nh(1)+nx(1):nh(end)+nx(end)]; title('Output of convolution y1(n)=x(n)*h(n)');
figure, xlim([ny(1)-1 ny(end)+1])
subplot(2,1,1) grid on
stem(ny,y);
CODE AND OUTPUT FOR TASK 3:
x=[0 0 1 0 0]; %let the input signal y3=conv(conv(x,h1),h2);

h1=[10 20 30 40 50] ; h2=[110 120 ny3=[nx(1)+nh1(1)+nh2(1):nx(end)+nh1(end)+n


h2(end)];
130 140 150];
figure,%distributive property
org_x=3;
title('Distributive property')
nx=(0:length(x)-1)-org_x+1 ;
subplot(2,1,1);
nh1=[0:length(h1)-1];
stem(ny,y);
nh2=(0:length(h2)-1);
xlim([ny(1)-1 ny(end)+1]);
%distributive property
xlabel('Time index n');
y=conv(x,h1+h2);
ylabel('Amplitude');
ny=[nx(1)+nh1(1):nx(end)+nh1(end)];
title('Output of convolution
y1=conv(x,h1) +conv(x,h2) ;
y(n)=x(n)*[h1(n)+h2(n)]')
ny1=[nx(1)+nh1(1):nx(end)+nh1(end)]; subplot(2,1,2);
%associative property
stem(ny1,y1);
y2=conv(x,conv(h1,h2)); xlim([ny1(1)-1 ny1(end)+1]);
ny2=[nx(1)+nh1(1)+nh2(1):nx(end)+nh1(end)+n xlabel('Time index n');
h2(end)];
ylabel('Amplitude'); title('Output of convolution
title('Output of convolution y(n)=x(n)*h1(n) + y(n)=x(n)*[h1(n)*h2(n)]');
x(n)*h2(n)'); subplot(2,1,2)
figure, %associative property stem(ny3,y3);
subplot(2,1,1); xlim([ny3(1)-1 ny3(end)+1]);
stem(ny2,y2); xlabel('Time index n');
xlim([ny2(1)-1 ny2(end)+1]); ylabel('Amplitude');
xlabel('Time index n'); title('Output of convolution
ylabel('Amplitude'); y(n)=[x(n)*h1(n)]*h2(n)');
OUTPUT FOR TASK 4:
OUTCOME BASED EDUCATION
LAB RUBRICS
Digital Signal Processing (EE-394)
Student Name: Roll No: Section:

Excellent Good Average Poor Comments


100% 75% 50% 0%
The background The concepts The theoretical There is no
knowledge of
Clarity of theory is fully are not knowledge or
the background
Concepts grasped. completely concepts are
theory/concepts.
cleared. inappropriate.
The logic of
MATLAB the
program/Simulink program/model
The The
model is properly is correct but
Obtaining program/model program/model
simulated to there is some
plot/results has some is not correct.
obtain the desired syntax error to
major errors.
response. obtain the
results

Lab task Late


Submission submission of
submitted in time.
Lab Task
LAB SESSION 04

OBJECTIVE:

To study discrete time correlation and apply it to real data to observe the correlation
between two signals.

THEORY:

1. Correlation is given as where ‘l’ is the lag. This is called cross-correlation and it
gives the magniyude and location of similarity between two signals. The
correlation between x(n) and y(n) . It is given as:

2. Generally rxy(l) = ryx(l). These two are the same when x(n) and y(n) are the same signals
or when x(n) and y(n) are even symmetric signals .

3. The length of the resulting correlation sequence is N+M-1, where N and M are the
lengths of the two signals.

4. Correlation may also be computed using convolution algorithm with a modification


that we need to fold one of the signals before applying convolution.
Mathematically, rxy(n)= x(n) * y(-n)

STEPS:

1. Generate two sinusoids of length 10 and fd = 0.1 with variable phase.


2. Apply correlation and check for certain properties such as magnitude and location
of maximum correlation with varying phases.

PROCEDURE:

1.Make a folder at desktop and name it as your current directory within MATLAB. -
2.Open M-file editor and write the following code: )
clear all;
close all;
clc;
n = [0:9];
ph1 = 0;
ph2 = 0;
x = sin(2*pi*0.1*n + ph1);
org_x = 1;
nx = [0 : length(x)-1]- org_x + 1;

y = sin(2*pi*0.1*n + ph2);
org_y = 1;
ny = [0 : length(y)-1]- org_y + 1;

rxy = xcorr(x,y);
nr = [nx(1)-ny(end) : nx(end)-ny(1)];

[maxR indR] = max(rxy);

disp(['The correlation at lag zero is: ' num2str(rxy(find(nr==0)))


'.']);
disp(['The maximum correlation is at lag ' num2str(nr(indR)) '.']);
|
figure,
subplot(3,1,1),
stem(nx,x);
xlabel('Time index n');
ylabel('Amplitude');
xlim([nx(1)-1 nx(end)+1]);
title('Signal x(n)');
grid;

subplot(3,1,2),
stem(ny,y);
xlabel('Time index n');
ylabel('Amplitude');
xlim([ny(1)-1 ny(end)+1]);
title('Signal y(n)');
grid;

subplot(3,1,3)
stem(nr,rxy);
xlabel('Time index n');
ylabel('Amplitude');
xlim([nr(1)-1 nr(end)+1]);
title('Cross Correlation');
grid;

3.Save the file as P041.m in your current directory and ‘run’ it.
Learn the specific logical bits of the code and make notes

Now modify the phase of the second signal to pi/2 (it will make it cosine) and observe
the correlation at lag zero. Modify the phase again to ‘pi’ and observe.

1. Check for auto-correlation (ph1 = ph2) that the lag zero value gives the energy of the
Signal.
2. Observe that the commutative property does not hold.

RESULT:
CODING
The correlation at lag zero is: 4.1633e-16.
The maximum correlation is at lag 2.
close all;
clc;
n = [0:9];
ph1 = 0;
ph2 = pi/2;
x = sin(2*pi*0.1*n + ph1);
org_x = 1;
nx = [0 : length(x)-1]- org_x + 1;
y = sin(2*pi*0.1*n + ph2);
org_y = 1;
ny = [0 : length(y)-1]- org_y + 1;
rxy = xcorr(x,y);
nr = [nx(1)-ny(end) : nx(end)-ny(1)];
[maxR indR] = max(rxy);
disp(['The correlation at lag zero is: ' num2str(rxy(find(nr==0))) '.']);
disp(['The maximum correlation is at lag ' num2str(nr(indR)) '.']);
figure,
subplot(3,1,1),
stem(nx,x);
xlabel('Time index n');
ylabel('Amplitude');
xlim([nx(1)-1 nx(end)+1]);
title('Signal x(n)');
grid;
subplot(3,1,2),
stem(ny,y);
xlabel('Time index n');
ylabel('Amplitude');
xlim([ny(1)-1 ny(end)+1]);
title('Signal y(n)');
grid;
subplot(3,1,3)
stem(nr,rxy);
xlabel('Time index n');
ylabel('Amplitude');
xlim([nr(1)-1 nr(end)+1]);
title('Cross Correlation');
grid;
Please write in exercise book.
EXERCISE:
1. Now modify the phase of the second signal to pi/2 (it will make it cosine)and observe the
correlation at lag zero.
2. Modify the phase again to ‘pi’ and observe.
3. Check for auto-correlation (ph1 = ph2) that the lag zero value gives the m energy of the
signal.
4. Observe that the commutative property does not hold.
5. Modify the code, such that the correlation is obtained using convolution command.
6. Calculate correlation between voltages of any two phases of a 10HP motor Using the data
given below. First use Ms. Excel to copy data and then calculate correlation.

CODE FOR THE TASK 1:


n=[0:9]; xlim([nx(1)-1 nx(end)+1]);
ph1=0; title('Signal x(n)');
ph2=pi/2; grid;
x=sin(2*pi*0.1*n + ph1) ; subplot(3,1,2),
y=sin(2*pi*0.1*n + ph2) ; stem(ny,y);
nx=[0:length(x)-1]; xlabel('Time index n');
ny=[0:length(y)-1]; ylabel('Amplitude');
rxy=xcorr(x,y); xlim([ny(1)-1 ny(end)+1]);
nr=[nx(1)-ny(end):nx(end)-ny(1)]; title('Signal y(n)');
[maxR indR]=max(rxy); grid;
disp(['The correlation at lag zero is:' subplot(3,1,3)
num2str(rxy(find(nr==0)))]); stem(nr,rxy);
figure, xlabel('Time index n');
subplot(3,1,1), ylabel('Amplitude');
stem(nx,x); xlim([nr(1)-1 nr(end)+1]);
xlabel('Time index n'); title('Cross Correlation');
ylabel('Amplitude'); grid;
The correlation at lag zero is:4.1633e-16
CODE FOR THE TASK 2:

n=[0:9]; xlim([nx(1)-1 nx(end)+1]);


ph1=0; title('Signal x(n)');
ph2=pi; grid;
x=sin(2*pi*0.1*n + ph1) ; subplot(3,1,2),
y=sin(2*pi*0.1*n + ph2) ; stem(ny,y);
nx=[0:length(x)-1]; xlabel('Time index n');
ny=[0:length(y)-1]; ylabel('Amplitude');
rxy=xcorr(x,y); xlim([ny(1)-1 ny(end)+1]);
nr=[nx(1)-ny(end):nx(end)-ny(1)]; title('Signal y(n)');
[maxR indR]=max(rxy); grid;
Ex=abs(rxy(find(nr==0))).^2; subplot(3,1,3)
disp(['The correlation at lag zero is:' stem(nr,rxy);
num2str(rxy(find(nr==0)))]); xlabel('Time index n');
figure, ylabel('Amplitude');
subplot(3,1,1), xlim([nr(1)-1 nr(end)+1]);
stem(nx,x); title('Cross Correlation');
xlabel('Time index n'); grid;
ylabel('Amplitude');
The correlation at lag zero is:-5

CODE FOR TASK 3:


n=[0:9]; ny=[0:length(y)-1];
ph1=pi/2; rxy=xcorr(x,y);
ph2=pi/2; nr=[nx(1)-ny(end):nx(end)-ny(1)];
x=sin(2*pi*0.1*n + ph1) ; disp(['The auto correlation gives '
y=sin(2*pi*0.1*n + ph2) ; num2str(Ex) 'J energy.']);
nx=[0:length(x)-1];

OUTPUT:
The auto correlation gives 25J energy.

CODE FOR TASK 4:


n=[0:9]; figure,
ph1=0; subplot(2,1,1);
ph2=3*pi/2; stem(nrxy,rxy);
x=sin(2*pi*0.1*n + ph1) ; xlabel('Time index n');
y=sin(2*pi*0.1*n + ph2) ; ylabel('Amplitude');
nx=[0:length(x)-1]; xlim([nrxy(1)-1 nrxy(end)+1]);
ny=[0:length(y)-1]; subplot(2,1,2);
rxy=xcorr(x,y); stem(nryx,ryx);
nrxy=[nx(1)-ny(end):nx(end)-ny(1)]; xlabel('Time index n');
ryx=xcorr(y,x); ylabel('Amplitude');
nryx=[nx(1)-ny(end):nx(end)-ny(1)]; xlim([nryx(1)-1 nryx(end)+1]);
After observing stems it can be clearly seen that the cross correlation does
not hold commutative property rather than after applying commutative
property we can clearly seen that the two plots in cross correlation becomes
folded version of each other.

CODE FOR TASK 5:


n=[0:9]; title('Signal x(n)');
ph1=0; xlim([nx(1)-1 nx(end)+1]);
ph2=0; subplot(3,1,2);
x=sin(2*pi*0.1*n + ph1); stem(ny,y);
y=sin(2*pi*0.1*n + ph2) ; xlabel('Time index n');
yf=fliplr(y) ; ylabel('Amplitude');
nx=(0:length(x)-1); title('Signal y(n)');
ny=(0:length(y)-1); xlim([ny(1)-1 ny(end)+1]);
nyf=(0:length(yf)-1); subplot(3,1,3)
rxy=conv(x,yf); stem(nr,rxy);
nr=[nx(1)-nyf(end):nx(end)-nyf(1)]; xlabel('Time index n');
figure, ylabel('Amplitude');
subplot(3,1,1); title('Correlation between x(n) and y(n)');
stem(nx,x); xlim([nr(1)-1 nr(end)+1]);
xlabel('Time index n');
ylabel('Amplitude');

Voltage A Min Voltage B Min Voltage C Min


189.358 153.917 195.735
189.175 159.719 201.877
188.783 161.575 186.718
188.757 172.186 187.659
176.995 173.206 205.876
180.472 176.865 204.831
180.524 176.917 192.494
180.262 189.28 199.839
181.778 189.828 211.887
179.975 189.462 211.94
178.642 189.253 212.462
180.315 188.94 193.749
180.707 190.377 200.492
180.262 190.194 201.433
180.628 190.064 202.635
180.315 189.907 200.701
179.635 189.541 203.289
179.243 189.567 202.635
179.4 189.619 200.989
180.576 189.044 197.591
180.837 189.123 199.865
180.184 189.332 201.093
180.08 189.097 201.041
177.675 189.044 199.656
175.297 189.018 198.558
173.99 189.123 204.595

CODE FOR TASK 6:


VA=[189.358 189.175 188.783 188.757 176.995 180.472 180.524 180.262 181.778 179.975
178.642 180.315 180.707 180.262 180.628 180.315 179.635 179.243 179.4 180.576 180.837
180.184 180.08 177.675 175.297 173.99 ];
VB=[153.917 159.719 161.575 172.186 173.206 176.865 176.917 189.28 189.828 189.462
189.253 188.94 190.377 190.194 190.064 189.907 189.541 189.567 189.619 189.044 189.123
189.332 189.097 189.044 189.018 189.123 ] ;
nVA=[0:length(VA)-1]; stem(nVB,VB);
nVB=[0:length(VB)-1]; xlabel('Time index n');
rxy=xcorr(VA,VB); ylabel('Voltage VB');
nr=[nVA(1)-nVB(end):nVA(end)-nVB(1)]; title('Signal VB');
figure, xlim([nVB(1)-1 nVB(end)+1]);
subplot(3,1,1); subplot(3,1,3);
stem(nVA,VA); stem(nr,rxy);
xlabel('Time index n'); xlabel('Time index n');
ylabel('Voltage VA'); ylabel('Amplitude');
title('Signal VA'); title('Correlation between VA and VB');
xlim([nVA(1)-1 nVA(end)+1]); xlim([nr(1)-1 nr(end)+1]);
subplot(3,1,2);
OUTCOME BASED EDUCATION
LAB RUBRICS
Digital Signal Processing (EE-394)
Student Name: Roll No: Section:

Excellent Good Average Poor Comments


100% 75% 50% 0%
The background The concepts The theoretical There is no
knowledge of
Clarity of theory is fully are not knowledge or
the background
Concepts grasped. completely concepts are
theory/concepts.
cleared. inappropriate.
The logic of
MATLAB the
program/Simulink program/model
The The
model is properly is correct but
Obtaining program/model program/model
simulated to there is some
plot/results has some is not correct.
obtain the desired syntax error to
major errors.
response. obtain the
results

Lab task Late


Submission submission of
submitted in time.
Lab Task
LAB SESSION 05

OBJECTIVE:
To study the computer implementation of Discrete Fourier T transform and Inverse Fourier
Transform using Twiddle factor.

THEORY:
The formulas for the DFT and IDFT are given as
X(K) = ∑ −1
( )
; k=0,1,……N-1
=0

X(n) = ; k=0,1,……N-1
1 ∑ −1 ( ) −

=0

Where by definition WN =
− 2

Which is an Nth root of unity. Where WN is called Twiddle Factor , also


[WN] = [ − 2 / ]kn

W = − 2 /

DFT analysis equation in matrix form is


XN=[ W ]xN

DFT synthesis equation in matrix form is


-1
xN= [W ] XN

PROCEDURE:

TASK
Compute 4 point DFT of x(n)= ( 1,2,3,0).

STEPS
1.Generate given sequence in Matlab .
2.Take N-=4 to calculate 4-point DFT.
3.Define 0: N-1 point vector for time and frequency samples.
4.Define W matrix and then use DFT analysis equation to compute DFT.

close all, k=[0:1:N-1];


clear all; WN=exp(-j*2*pi/N);
clc; nk=n'*k;
x=[1 ,2 ,3 ,0]; WNnk=WN.^nk;
N=4; Xk=x*WNnk
n=[0:1:N-1];
LAB TASK

Task-1: Forward DFT using matrices


Develop a MATLAB code to find the Forward DFT output of the following time domain sequence by using DFT equation in matrix form. Also plot the magnitude and phase spectrum. Take =
1000 ⁄

( )={ . , 0.3535, 0.6464, 1.0607 , 0.3535, − 1.0607, − 1.3535 , − 0.3535}

Task-2: Inverse DFT using Matrix inversion


Develop a MATLAB code to find the inverse DFT output of the following frequency domain sequence by
using iDFT equation in matrix form (use matrix inversion).

0, 0,
( ) = { , 4∠ − 90 ∘, 2∠45 ∘, 0, 2∠−45∘, 4∠90 ∘}

Task-3: Inverse DFT using Conjugate method


Develop a MATLAB code to find the inverse DFT output of the following frequency domain sequence by
using iDFT equation in matrix form (use conjugate method).

0, 0,
( ) = { , 4∠ − 90 ∘, 2∠45 ∘, 0, 2∠−45∘, 4∠90 ∘}

CODE FOR TASK:


close all, a=4*exp(1i*deg2rad(-90));
clear all; b=2*exp(1i*deg2rad(45));
clc; c=exp(1i*deg2rad(-45));
d=4*exp(1i*deg2rad(45)); nk=n'*k;
x=[0,a,b,0,0,0,c,d];
N=8; WNnk=WN.^nk;
n=[0:1:N-1];
k=[0:1:N-1]; Xk=x*WNnk
f=conj(WNnk)
WN=exp(-j*2*pi/N); xk=(Xk*f)/N
OUTCOME BASED EDUCATION
LAB RUBRICS
Digital Signal Processing (EE-394)
Student Name: Roll No: Section:

Excellent Good Average Poor Comments


100% 75% 50% 0%
The background The concepts The theoretical There is no
knowledge of
Clarity of theory is fully are not knowledge or
the background
Concepts grasped. completely concepts are
theory/concepts.
cleared. inappropriate.
The logic of
MATLAB the
program/Simulink program/model
The The
model is properly is correct but
Obtaining program/model program/model
simulated to there is some
plot/results has some is not correct.
obtain the desired syntax error to
major errors.
response. obtain the
results

Lab task Late


Submission submission of
submitted in time.
Lab Task
LAB SESSION 06

OBJECTIVE:
To observe/find different frequency components in an audio signal and plot it with different x_
axes .

THEORY:

PROCEDURE:

1. Load an audio file ‘noisy.wav’ into Matlab.


2. There is a tone added to the speech in this file. The objective is to find the frequency of
this tone.
3. Computing the DFT of this signal;
4. Generating frequency vector in Hz.
5. Displaying the DFT and observing the frequency of the added tone.

STEPS
1. Make a folder at desktop and name it as your current directory within MATLAB.
2. Copy the audio file ‘noisy.wav’ into your current directory.
3. Open M file editor and write the following code:

clear all; clc; close all;


[y,Fs] = wavread('noisy.wav');
Ts = 1/Fs;
n = [0:length(y)-1];
t = n.*Ts; k = n; Df
= Fs./length(y); F =
k.*Df;
Y = fft(y);
magY = abs(Y);
sound(y,Fs);
figure,

subplot(2,1,1);
plot(F,magY);
grid on;
xlim([0 Fs/2]);
xlabel('Frequency (Hz)');
ylabel('DFT Magnitude');
title('Discrete Fourier Transform');

subplot(2,1,2);
plot(F,magY);
grid on;

xlim([0 2000]);
xlabel('Frequency (Hz)');
ylabel('DFT Magnitude');
title('Discrete Fourier
Transform');

4. Save the file as P081.m in your current directory and run it.

RESULT:

Explore and take notes.


EXERCISE:

Use recorded data,


1. Plot different frequencies present in it with
a) x-axis as time
b )x-axis as frequency. (Take FFT and plot).
2. Calculate the amount of energy present in fundamental frequency.
3. Calculate the amount of energy present in different harmonics.

CODE FOR THE TASK:


[x,Fs]=audioread('output.wav'); num2str(i-1) ' is: ' num2str(E_(i))]);
sound(x,Fs); end
n=[0:length(x)-1] ; i=i+1;
t=n.*1/Fs; end
k=[0:length(x)-1]; figure,
DF=Fs/length(x) ; subplot(2,1,1);
F=k.*DF ; plot(F,mag) ;
X=fft(x); xlabel('Frequency in Hz');
mag=abs(X) ; ylabel('magnitude');
i=1; xlim([0 2000]);
while i<=5 grid on
E_(i)=[x(i)].^2; subplot(2,1,2);
if i==1 plot(t,x) ;
disp(['The amount of energy at fundamental xlabel('Time domain (sec)');
harmonic is: ' num2str(E_(i))]); ylabel('magnitude');
else xlim([0 6]);
disp(['The amount of energy at harmonic ' grid on
OUTCOME BASED EDUCATION
LAB RUBRICS
Digital Signal Processing (EE-394)
Student Name: Roll No: Section:

Excellent Good Average Poor Comments


100% 75% 50% 0%
The background The concepts The theoretical There is no
knowledge of
Clarity of theory is fully are not knowledge or
the background
Concepts grasped. completely concepts are
theory/concepts.
cleared. inappropriate.
The logic of
MATLAB the
program/Simulink program/model
The The
model is properly is correct but
Obtaining program/model program/model
simulated to there is some
plot/results has some is not correct.
obtain the desired syntax error to
major errors.
response. obtain the
results

Lab task Late


Submission submission of
submitted in time.
Lab Task
LAB SESSION 07

OBJECTIVE:
The purpose of this lab is to become familiar with the practical constraint in calculating the Fourier
transform of any real-world signal. i.e. DFT or FFT leakage and its solution.
INTRODUCTION:

When DFT Leakage will happen?


If the input signal contains frequency which is integer multiple of the DFT frequency
resolution (Fs/N) or analysis frequencies (k*Fs/N), it will show up only at the correct
DFT output bin or frequency. All other DFT output bins will be zero. It means
spectrum is correctly estimated/calculated and no DFT leakage occurs in this case. But
this is an ideal case to dictate the ingredient components of real-world signal to become
integral multiple of DFT frequency resolution or analysis frequencies.
If the input signal contains any frequency which is NOT integer multiple of analysis
frequencies (which is the most probably happening case), it will show up in all DFT
output bins or leaks out to all DFT output bins. This effect is simply known as DFT
Leakage/ Spectral Leakage. We typically say that input signal energy shows up in all
of the DFT's output bins. Think of bins as buckets to which you pour energy from
infinite amount of possible frequency components. If the signal's frequency is not
matched with exact frequency bin then it smears leaks over rest of the bins.

DFT Leakage and Use of Window functions:

Another way to think about DFT leakage is, when the input signal has frequencies which do
not complete integer number of cycles in the sample interval N (Length of DFT), the input to
DFT seems like abruptly starting and abruptly ending giving rise to side lobes in the DFT
output. A selection of finite input samples means that we have multiplied the input sequence
by a rectangular window. The DFT output will follow the Sinc shape which is the Fourier
Transform of the rectangular function. For real life signals, DFT leakage can never be entirely
eliminated because no matter what point DFT you take and no matter what your sampling rate
is, there is a very high probability that the input would contain a frequency which is not an
integer multiple of analysis frequencies and DFT leakage would happen. We can minimize
DFT leakage by multiplying the input with a smooth window which will make the input go
from zero to maximum and back to zero in a very smooth and slow fashion. This reduces the
abruptness of the input which minimizes the higher sides lobes in the DFT output.
Some Window functions and their spectrums

LAB TASKS

Task-1: Condition of DFT/FFT Leakage

Fill the following table specifying whether DFT leakage would happen or not?

Sr. Fs N Input contains these frequencies in DFT Leakage


Hz
No. [samples/sec] [samples] Yes or No?
1 8000 8 1000
2 8000 9 1000, 2000
3 16000 32 250, 500, 1000, 2250
4 22050 128 11025
5 44100 4096 1000, 2000, 2500, 11025
6 44100 44100 1101, 2202, 3303
7 96000 48000 1200, 1202, 2002, 2003

Task-2: Calculation of cycles in N sample interval (length of DFT/FFT)

Calculate exactly how many cycles of the 440 Hz wave sampled at 1000 Hz are contained in the 1024-
sample window. Explain why this number indicates leads to spectral leakage.

Task-3: DFT or FFT Leakage

a) Generate two sinusoids of frequencies 500 Hz and 600 Hz. Both signals should be
sampled at 8000 samples/sec.
b) Take 64-point DFT of both signals and observe their spectrums.
c) Comment on the DFT Leakage for both cases.
PROCEDURE:

clear all; close all; clc; subplot(4,1,2)


stem(Fk,X1,'r','filled')
F1 = 500; xlabel('Frequency (Hz)')
F2 = 600; ylabel('DFT Magnitude')
Fs = 8000; title('Spectrum')
Ts = 1/Fs; xlim([0 Fs])
N=64; grid
n =[0:N-1];
k = n; subplot(4,1,3)
Df = Fs/N; plot(n*Ts,x2,'r')
Fk = k.*Df; hold
stem(n*Ts,x2,'filled' )
x1 = sin(2*pi*F1*n*Ts); xlabel('Time (sec)')
X1 = abs(fft(x1,N)); ylabel('Amplitude')
x2 = sin(2*pi*F2*n*Ts); title('Signal')
X2 = abs(fft(x2,N)); grid
axis tight
subplot(4,1,1)
plot(n*Ts,x1,'r') subplot(4,1,4)
hold stem(Fk,X2,'r','filled')
stem(n*Ts,x1,'filled' ) xlabel('Frequency (Hz)')
xlabel('Time (sec)') ylabel('DFT Magnitude')
ylabel('Amplitude') title('Spectrum')
title('Signal') xlim([0 Fs])
grid grid
axis tight
Task-4: Evaluation of Various Window functions using wintool of MATLAB or Window
Analysis & Designing

Type wintool (Window Design & Analysis Tool) in the command window of MATLAB. Observe
different window functions like Rectangular (no window), Triangular, Hamming, Hanning,
Blackman, Flat Top, Gaussian etc. Attach the time domain and frequency domain representations of
various window functions and make a comparison between their performances on the basis of following
attributes:

a)
b)
Leakage Factor in %
Width of the main lobe in fraction of π ( ⁄ ) or ⁄2 ( ⁄sec )

c) The first (highest)side lobe magnitude in dB


d) The rate of side lobes falloff in dB/octave

Task-5: Spectrum of Windowed signal


a) Generate a sinusoid of frequency 1050 Hz sampled at 8000 samples/sec. Initial phase is
⁄4.

b) Take 64-point DFT of the generated signal and observe the spectrum.
c) Comment on the DFT Leakage phenomenon for the above case.
d) Generate a triangular window function of 64 points.
e) Now apply a triangular window to the signal and then observe the spectrum of windowed
signal.
f) Change the window function from Triangular to Hamming, Hanning, Blackman, Flat
Top and Gaussian. Compare the results of these window functions on a signal's
spectrum.

PROCEDURE:

clear all;close all;clc;


subplot(3.2.3)
F = 1050; stem(n,w,'r','filled')
Fs = 8000; xlabel('Samples')
Ts = 1/Fs; ylabel('Amplitude')
N=64; title('Window Function')
n = [0:N-1]; grid
k = n; axis tight
Df = Fs/N;
Fk = k.*Df; subplot(3,2,4)
stem(abs(fft(w)),'r','filled')
x = sin(2*pi*F*n*Ts+pi/4); xlabel('Samples')
ylabel('Magnitude')
w = window(@triang,N); title('Spectrum of Window Function')
% w = hamming(N); grid
% w = hann(N); axis tight
% w = blackman(N);
% w= flattopwin(N); subplot(3,2,5)
% w= gausswin(N); plot(n*Ts,xw,'r')
hold
stem(n*Ts,xw,'filled' )
xw = x.*w'; % Windowed signal xlabel('Time (sec)')
X = abs(fft(x)); ylabel('Amplitude')
Xw = abs(fft(xw)); title('Windowed Signal')
grid
subplot(3,2,1) axis tight
plot(n*Ts,x,'r') subplot(3,2,6)
hold
stem(n*Ts,x,'filled') stem(Fk,Xw,'filled')
xlabel('Time (sec)') xlabel('Frequency (Hz)')
ylabel('Amplitude') ylabel('Windowed DFT Magnitude')
title('Signal') title('Windowed Spectrum')
grid grid
subplot(3,2,2)

stem(Fk,X,'filled')
xlabel('Frequency (Hz)')
ylabel('DFT Magnitude')
title('Spectrum')
grid

Triangle window:
Hamming Window:

Hanning Window:
Flat top Window:

Blackman Window:
Gaussian Window:

OUTCOME BASED EDUCATION


LAB RUBRICS
Digital Signal Processing (EE-394)
Student Name: Roll No: Section:

Excellent Good Average Poor Comments


100% 75% 50% 0%
The background The concepts The theoretical There is no
knowledge of
Clarity of theory is fully are not knowledge or
the background
Concepts grasped. completely concepts are
theory/concepts.
cleared. inappropriate.
The logic of
MATLAB the
program/Simulink program/model
The The
model is properly is correct but
Obtaining program/model program/model
simulated to there is some
plot/results has some is not correct.
obtain the desired syntax error to
major errors.
response. obtain the
results

Lab task Late


Submission submission of
submitted in time.
Lab Task
LAB SESSION 08

OBJECTIVE:
To study s-plane and plot impulse and frequency response for different pole zero location in s-
plane. Also to determine weather system is FIR or IIR.

THEORY:
The Laplace Transform of a general continuous time signal x (t) is defined as;

X(S) = ∫ x(t) e-st dt.


Where the complex variable s=δ+ j w, with δ and w the real and imaginary parts. CTFT is a
subset of Laplace when δ =0. Since ‘δ’ information is not present in CTFT, therefore information
about stability can only be obtained from Laplace. If pole lies on L.H.S of s-plane, system is
stable. If pole lies on R.H.S of s-plane, system is unstable. If pole lies on y(jw)-axis, system is
marginally stable or oscillatory. If system has FIR, it is stable. If system is IIR, it can be stable or
unstable.

PROCEDURE:
Generate pole zero constellation in s plane.
1. Plot corresponding Frequency (Bode magnitude) response.
2. Plot impulse response and determine that the system is FIR or IIR.
3. Modify location of poles in s plane to observe the corresponding change in frequency and
impulse response.

STEPS.
1. Make a folder at desktop and name it as your current directory within MATLAB.
2. Open M-file editor and write the following code:

clear all; pzmap(sys);


close all; xlim([-2 2]);
clc; ylim([-4 4]);
Num = poly([(0- subplot(3,1,2);
(i*(pi/2))),(0+(i*(pi/2)) [mag phase w]=bode(sys);
)]); mag=squeeze(mag);
Zeros=roots(Num) plot(w,mag);
Den = poly([-1,-1]); subplot(3,1,3);
poles=roots(Den) impulse(sys);
sys=tf(Num,Den) H=dfilt.df1(Num,Den);A=is
figure; fir(H)
subplot(3,1,1);
3. Save the file as P091.m in your current directory and ‘run’ it.

RESULT:

1. Learn the specific logical bits of the code and make notes.
2. Observe the plots.
3. Now, explain (write) in your own words the cause and effects of what you just saw.
EXERCISE:

Task-1: Analysis of 1st order Analog system [one pole at origin]

a) Generate pole zero constellation of an analog system in s-plane having only one pole at
s=0
b) Write the transfer function of a system
c) Plot the corresponding frequency response and impulse response of a system. Also
comment on the stability of system

PROCEDURE:

a) Make a folder at desktop and name it as your current directory within MATLAB.
b) Open M-file editor and write the following code:

clear all; close all; clc; subplot(3,1,2)


Num = [1]; [H w]=freqs(Num,Den);
Den = poly([0]); % At origin mag=abs(H);
plot(w,mag)
sys=tf(Num,Den) xlabel('Frequency (\omega)')
zeros=roots(Num) ylabel('Gain |H(\omega)|' )
poles=roots(Den) title('Frequency Response')
axis tight
figure; grid
subplot(3,1,1) subplot(3,1,3)
pzmap(sys)
xlim([-3 3]) impulse(sys,'r')
ylim([-3 3]) grid
Task-2: Analysis of 1st order Analog system [one pole at LHP]

a) Generate pole zero constellation of an analog system in s-plane having only one pole at: s = -0.5
b) Write the transfer function of a system
c) Plot the corresponding frequency response and impulse response of a system. Also comment on
the stability of system.

Task-3: Analysis of 1st order Analog system [one pole at RHP]

a) Generate pole zero constellation of an analog system in s-plane having only one pole at: s =0.5
b) Write the transfer function of a system
c) Plot the corresponding frequency response and impulse response of a system. Also comment on
the stability of system.

Task-4: Analysis of 2nd order Analog system [poles at j ]

a) Generate pole zero constellation of an analog system in s-plane having pure imaginary poles at: s= = ± ⁄2
b) Write the transfer function of a system
c) Plot the corresponding frequency response and impulse response of a system. Also comment on
the stability of system.

Task-5: Analysis of 2nd order Analog system [complex poles at LHP]


a) Generate pole zero constellation of an analog system in s-plane having complex conjugate poles at: = −0.5 ± ⁄2
b) Write the transfer function of a system
c) Plot the corresponding frequency response and impulse response of a system Also comment on
the stability of system.

Task-6: Analysis of 2nd order Analog system [complex poles at RHP]


a) Generate pole zero constellation of an analog system in s-plane having complex conjugate poles at: = +0.5 ± ⁄2
b) Write the transfer function of a system
c) Plot the corresponding frequency response and impulse response of a system. Also comment on
the stability of system.

Task-7: Analysis of 2nd order Analog system [complex zeros and poles at LHP]
a) Generate pole zero constellation of an analog system in s-plane for given roots. [zeros: = ± ( /2) & poles: = −0.2 ± ( /4)]
b) Write the transfer function of a system
c) Plot the corresponding frequency response and impulse response of a system. Also comment on
the stability of system.

Task-8: Effect of system's poles on system stability

Change the location of poles of a system defined in Task-7 from L.H.S of s-plane to axis first, and
then to R.H.S of s-plane and observe the effects on impulse response and frequency response of a
system.

Task-9: Effect of system's zeros on system stability


Modify the location of zeros of a system defined in Task-7 from = ± ( /2) to = −0.1 ± ( /2) , and then to = 0.1 ± ( /2) . Do not change the
location of poles. Does the impulse response change? Can we place analog system's zeros on the right hand plane [RHP]?

OUTCOME BASED EDUCATION


LAB RUBRICS
Digital Signal Processing (EE-394)
Student Name: Roll No: Section:

Excellent Good Average Poor Comments


100% 75% 50% 0%
The background The concepts The theoretical There is no
Clarity of
knowledge of
theory is fully are not knowledge or
the background
Concepts grasped. completely concepts are
theory/concepts.
cleared. inappropriate.
The logic of
MATLAB the
program/Simulink program/model
The The
model is properly is correct but
Obtaining program/model program/model
simulated to there is some
plot/results has some is not correct.
obtain the desired syntax error to
major errors.
response. obtain the
results

Lab task Late


Submission submission of
submitted in time.
Lab Task
LAB SESSION 09
OBJECTIVE:

To study z-plane and plot impulse and frequency response for different pole zero location in z-
plane.Also to determine weather system is FIR or IIR.

THEORY:

The z - Transform of a general discrete time



signal x(n) is defined as;
X (z) = ∑ = ( ) z-n

Where the complex variable z=r ∠w , with r the radius and w the angle. DTFT is a subset of z transform
when r =1. Since ‘r’ information is not present in DTFT, therefore information about stability in discrete
time can only be obtained from z transform. If pole lies inside the unit circle, system is stable. If pole lies
outside the unit circle, system is unstable. If pole lies at the unit circle, system is marginally stable or
oscillatory. If system has FIR, it is stable. If system is IIR, it can be stable or unstable .

PROCEDURE:

1. Generate pole zero constellation in z plane.


2. Plot corresponding Frequency (Bode magnitude) response.
3. Plot impulse response and determine that the system is FIR or IIR.
4. Modify location of poles in z plane to observe the corresponding change in frequency
and impulse response.

STEPS:

1. Make a folder at desktop and name it as your current directory within MATLAB.
2. Open M-file editor and write the following code:

clear all; Num1 = poly([j,-j]);


close all; Den1 = poly([exp(-
clc; 1),exp(-1)]);
Num = poly([(0- sys1=tf(Num1,Den1,1)
(i*(pi/2))),(0+(i*(pi/2))
)]); figure;
subplot(3,1,1);
Den = poly([-1,-1]); pzmap(sys1);
xlim([-2 2]);
ylim([-4 4]); subplot(3,1,3);
subplot(3,1,2); impulse(sys1);
[mag phase w]=bode(sys1); H=dfilt.df1(Num,Den);
mag=squeeze(mag); A=isfir(H)
figure;
plot(w,mag); pzmap(sys1)
xlim([0 100]) grid on;
3. Save the file as P010.m in your current directory and ‘run’ it.

RESULT:

1 Learn the specific logical bits of the code and make notes.
2 Observe the plots.
3 Now, explain (write) in your own words the cause and effects of what you just saw.
EXERCISE:
Task-1: Analysis of 1st order Digital system [one pole inside the Unit circle at DC]
a) Generate pole zero constellation in z-plane for the following roots of a digital system: [zero: = 0 & pole: = 0.9∠0]
b) Write the transfer function of a system [see command window]
c) Determine that whether the digital system is FIR or IIR
d) Plot corresponding frequency response and impulse response of a digital system

PROCEDURE:

c) Make a folder at desktop and name it as your current directory within MATLAB.
d) Open M-file editor and write the following code:

clear all;close all;clc; mag=abs(H);


plot((w*180)/pi,mag)
Num = poly([0]); xlabel('\omega (in degrees)')
Den = poly([0.9*exp(j*0)]); ylabel('Gain |H(\omega)|')
title('Frequency Response')
sys=tf(Num,Den,1) xlim([0 180])
%tf(num,den,Ts) for DT system where grid
Ts is sampling time subplot(3,1,3)
figure impulse(sys,'r')
subplot(3,1,1) grid
pzmap(sys)
xlim([-3 3]), ylim([-1.5 1.5]) H=dfilt.df1(Num,Den) %Digital Filter
subplot(3,1,2) with Direct Form-I approach
[H w]=freqz(Num,Den); A=isfir(H)
Task-2: Analysis of 1st order Digital system [pole inside the Unit circle at Fs/2]
a) Generate pole zero constellation in z-plane for the following roots of a digital system: [zero: = 0 & pole: = 0.9∠ ± ]
b) Write the transfer function of a system
c) Determine that whether the digital system is FIR or IIR
d) Plot corresponding frequency response and impulse response of a digital system

Task-3: Analysis of 1st order Digital system [one pole at the DC location of Unit circle]
a) Generate pole zero constellation in z-plane for the following roots of a digital system: [zero: = 0 & pole: = 1∠0]
b) Write the transfer function of a system
c) Determine that whether the digital system is FIR or IIR
d) Plot corresponding frequency response and impulse response of a digital system

Task-4: Analysis of 1st order Digital system [pole at the Fs/2 location of Unit circle]
a) Generate pole zero constellation in z-plane for the following roots of a digital system: [zero: = 0 & pole: = 1∠ ± ]
b) Write the transfer function of a system
c) Determine that whether the digital system is FIR or IIR
d) Plot corresponding frequency response and impulse response of digital system

Task-5: Analysis of 1st order Digital system [one pole outside the Unit circle at DC]
a) Generate pole zero constellation in z-plane for the following roots of a digital system: [zero: = 0 & pole: = 1.1∠0]
b) Write the transfer function of a system
c) Determine that whether the digital system is FIR or IIR
d) Plot corresponding frequency response and impulse response of a digital system

Task-6: Analysis of 1st order Digital system [pole outside the Unit circle at Fs/2]
a) Generate pole zero constellation in z-plane for the following roots of a digital system: [zero: = 0 & pole: = 1.1∠ ± ]
b) Write the transfer function of a system
c) Determine that whether the digital system is FIR or IIR
d) Plot corresponding frequency response and impulse response of a digital system
Task-7: Analysis of 2nd order Digital system [poles at origin]
a) Generate pole zero constellation in z-plane for the following roots of a digital system: [zeros: = 0.8944∠ ± 2 /3) & poles: = 0]
b) Write the transfer function of a system
c) Determine that whether the digital system is FIR or IIR
d) Plot corresponding frequency response and impulse response of a digital system

Task-8: Analysis of 2nd order Digital system [complex poles inside the Unit circle]
a) Generate pole zero constellation in z-plane for the following roots of a digital system: [zeros: = 1∠ ± /2) & poles: = 0.8∠ ± /4]
b) Write the transfer function of a system
c) Determine that whether the digital system is FIR or IIR
d) Plot corresponding frequency response and impulse response of a digital system

Task-9: Analysis of 2nd order Digital system [complex poles at the Unit circle]
a) Generate pole zero constellation in z-plane for the following roots of a digital system: [zeros: = 1∠ ± /2) & poles: = 1∠ ± /4]
b) Write the transfer function of a system
c) Determine that whether the digital system is FIR or IIR
d) Plot corresponding frequency response and impulse response of a digital system

Task-10: Analysis of 2nd order Digital system [complex poles outside the Unit circle]
a) Generate pole zero constellation in z-plane for the following roots of a digital system: [zeros: = 1∠ ± /2) & poles: = 1.2∠ ± /4]
b) Write the transfer function of a system
c) Determine that whether the digital system is FIR or IIR
d) Plot corresponding frequency response and impulse response of a digital system

Task-11: Analysis of 2nd order Digital system [complex zeros inside the Unit circle]
a) Generate pole zero constellation in z-plane for the following roots of a digital system: [zeros: = 0.8∠ ± /2) & poles: = 0.8∠ ± /4]
b) Write the transfer function of a system
c) Determine that whether the digital system is FIR or IIR
d) Plot corresponding frequency response and impulse response of a digital system

Task-12: Analysis of 2nd order Digital system [complex zeros outside the Unit circle]
a) Generate pole zero constellation in z-plane for the following roots of a digital system: [zeros: = 1.2∠ ± /2) & poles: = 0.8∠ ± /4]
b) Write the transfer function of a system
c) Determine that whether the digital system is FIR or IIR
d) Plot corresponding frequency response and impulse response of a digital system.
OUTCOME BASED EDUCATION
LAB RUBRICS
Digital Signal Processing (EE-394)
Student Name: Roll No: Section:

Excellent Good Average Poor Comments


100% 75% 50% 0%
The background The concepts The theoretical There is no
knowledge of
Clarity of theory is fully are not knowledge or
the background
Concepts grasped. completely concepts are
theory/concepts.
cleared. inappropriate.
The logic of
MATLAB the
program/Simulink program/model
The The
model is properly is correct but
Obtaining program/model program/model
simulated to there is some
plot/results has some is not correct.
obtain the desired syntax error to
major errors.
response. obtain the
results

Lab task Late


Submission submission of
submitted in time.
Lab Task
LAB SESSION 10(a)
OBJECTIVE:
Object of this lab is introduction to digital filters and its types, design FIR filter and study how it
performs filtering on a signal. Further truncate different types of FIR filter like Low Pass, High
Pass, Band Pass using different windows like rectangular, Kaiser Etc. and compare the results
obtained from different windows.

THEORY:
The process of deriving a realizable transfer function of a digital filter by considering given
frequency response specifications is known as digital filter design. The digital filter can be
classified as:

1. Finite –duration impulse response (FIR) filter


2. Infinite –duration impulse response (IIR) filter
In MATLAB, there are built in functions which can be used to design digital filter like IIR and
FIR.

The different types of FIR filters are listed as follows:

Window techniques based FIR filter.


1. Rectangular windows.
2. Hamming window.
3. Hanning window.
4. Blackman window.
5. Barlett window.
6. Kaiser window.
Equiripple linear phase FIR filter.
Least square error FIR filters.

The different types of IIR filters are listed as follows:

Butterworth filter
Chebyshev Type I filter
Chebyshev Type II
filter Elliptic filter

FIR digital filter operates on digital sample values. It uses current and past input samples to
produce a current output sample. It does not use previous output samples. There are various
types of FIR filter based on need viz. low pass, high pass, band pass and band stop, Low pass
filter.
Following points are usually considered to design FIR filter other the window type.
INPUT:
• Window Type
• Passband and stopband ripples
• passband and stopband edge frequencies
• sampling frequency
• order of the filter
• window coefficients

OUTPUT:
• magnitude and phase responses

COMPARISON OF FIR AND IIR FILTERS


1. FIR filters are Finite Impulse Response filters with no feedback, whereas IIR contains
feedback.
2. Transfer function of FIR filter does not contain any non-trivial poles. Their frequency
response is solely dependent on zero locations. IIR filters contain poles as well as zeros.
3. As there are no poles, FIR filters cannot become unstable; however, IIR filters can
become unstable if any pole lies outside the unit circle in z-plane.
4. More number of coefficients is needed to design the desired filter in FIR than IIR.

PROCEDURE:
TASK-1

1. Create a signal vector containing two frequencies as:


i) 100 Hz. and ii) 150 Hz. with Fs = 1000 Hz.
2. Design two band pass FIR filters with 64 coefficients and with pass bands as i) 125 to
175 Hz. and ii) 75 to 125 Hz.
3. Use both filters on the created signal and observe their outputs.
4. Plot frequency responses and pole-zero constellations of both filters and note
observations.

close all; clear all; clc; t = [0 : 1/Fs : 1]; % Time


% Frequencies in Hz. Vector
F = Fs*[0:length(t)-
1]/length(t); % Frequency
Vector
F1 = 100; F2 = 150;
x =
% Sampling Frequency exp(j*2*pi*F1*t)+2*exp(j*2*pi
in samples/sec. *F2*t); % Signal Vector

Fs = 1000; bh = fir1( 64 , [125


175]/500); % filter coeffs. subplot(5,1,4),
bl = fir1( 64 , [75 plot(F,abs(hl)); xlim([0
125]/500); % filter Fs/2]); title('Frequency
coeffs. response of Filter Two');
[hh,wh]=freqz(bh,1,leng
th(t),'whole'); % subplot(5,1,5),
Frequency response for plot(F,abs(fft(yl)));
filter 1 xlim([0 Fs/2]);
[hl,wl]=freqz(bl,1,leng title('FFT of filtered
th(t),'whole'); % signal from filter two');
Frequency response for xlabel('Hz.')
filter 2
% Filter operation - see
filtfilt in help to learn
what it does % Pole Zero
Constellations
yh = filtfilt(bh,1,x);
yl = filtfilt(bl,1,x); [bh,ah] =
eqtflength(bh,1);
% Plotting [zh,ph,kh] =
tf2zp(bh,ah);
figure, subplot(5,1,1), [bl,al] =
eqtflength(bl,1);
plot(F,abs(fft(x))); [zl,pl,kl] =
tf2zp(bl,al);
xlim([0 Fs/2]);
title('FFT of figure,
original signal'); subplot(1,2,1),
pzplane(bh,ah);
subplot(5,1,2), xlim([-1.5 1.5]);
plot(F,abs(hh)); xlim([0 ylim([-1.5 1.5]);
Fs/2]); title('Frequency title('Filter_One');
response of Filter One'); subplot(1,2,2),
pzplane(bl,al);
subplot(5,1,3), xlim([-1.5 1.5]);
plot(F,abs(fft(yh))); ylim([-1.5 1.5]);
xlim([0 Fs/2]); title('Filter Two');
title('FFT of filtered
signal from filter one');

TASK -2
Write a program to design a FIR filter using Hanning windows,take inputs from user for design
values of filter.

close all;
clear all;
clc;

fp=input('Enter the pass band frequency');


fs=input('Enter the stop band frequency');
rp=input(' Enter the pass band attenuation');
rs=input('Enter the stop band attenuation');
f=input(' Enter the sampling frequency');

% Calculating filter order

num=-20*log10(sqrt(rp*rs))-13;
dem=14.6*(fs-fp)/f;
n=ceil(num/dem);
n=abs(n);

% Normalizing the frequencies

wp=2*fp/f;
ws=2*fs/f;
wn=(ws+wp)/2;

%Adjusting the filter order. The order of window must be an odd


number
%and the order of filter must be one less than that of the
window

if (rem(n,2)==0)
m=n+1;
else
m=n;
n=n-1;
end

%Window sequence calculation


w=hann(m);

%Calculation of filter coefficients


b=fir1(n,wn,'low',w);

%Plotting the filter response


freqz(b,1,500,3000);
TITLE('Magnitude and Phase response');

TASK-3
Write a program for FIR(Finite Impulse Response) filter like Low pass FIR filter, High pass FIR
filter, Band pass FIR filter and Band stop FIR filter using Rectangular window using MATLAB .
ALGORITHM:
LOW PASS FILTER:
Step 1: Read the input sequence
Step 2: Perform low pass filter calculations
Step 3: Plot the output sequences

HIGH PASS FILTER:


Step 1: Read the input sequence
Step 2: Perform high pass filter calculations
Step 3: Plot the output sequences

BAND PASS FILTER:


Step 1: Read the input sequence
Step 2: Perform band pass filter calculations
Step 3: Plot the output sequences

BAND STOP FILTER:


Step 1: Read the input sequence
Step 2: Perform band stop filter calculations
Step 3: Plot the output sequences

PROGRAM:
clc; if(rem(n,2)~=0)
clear all; n1=n;
close all; n=n-1;
rp=input('Enter the passband end
ripple(rp):'); y=boxcar(n1);
rs=input('Enter the stopband %Low pass filter
ripple(rs):'); b=fir1(n,wp,y);
[h,o]=freqz(b,1,256);
fp=input('Enter the passband
m=20*log10(abs(h));
frequency(fp):'); subplot(2,2,1);
fs=input('Enter the stopband plot(m);
frequency(fs):'); ylabel('Gain(db)->');
f=input('Enter the sampling xlabel('(a)Normalised
frequency(f):'); frequency->');
wp=2*fp/f; %High pass filter
ws=2*fs/f; b=fir1(n,wp,'high',y);
num=-20*log10(sqrt(rp*rs))- [h,o]=freqz(b,1,256);
13; m=20*log10(abs(h));
dem=14.6*(fs-fp)/f;
n=ceil(num/dem); subplot(2,2,2);
plot(m);
n1=n+1;
ylabel('Gain(db)');
xlabel('(b)Normalised %Band stop
frequency'); filter==============
%Band pass filter wn=[wp*ws];
wn=[wp*ws]; b=fir1(n,wn,'stop',y);
b=fir1(n,wn,y); [h,o]=freqz(b,1,256);
[h,o]=freqz(b,1,256); m=20*log10(abs(h));
m=20*log10(abs(h)); subplot(2,2,4);
subplot(2,2,3); plot(m);
plot(m); ylabel('Gain(db)');
ylabel('Gain(db)'); xlabel('(d)Normalised
xlabel('(c)Normalised frequency-');
frequency');

EXERCISE:
Task-1: Designing of an IIR Butterworth filter

Design a 16th order IIR Butterworth bandpass filter with the following specifications:
a) Normalized pass band edges at 0.40 and 0.65
b) Normalized stop band edges at 0.3 and 0.75
c) Pass band ripple 1 dB
d) Minimum stop band attenuation 40 dB
Attach (i) Magnitude response (ii) Phase response (iii) Impulse Response (iv) Pole Zero plot

Task-2: Designing of an IIR Chebyshev Type II filter


Design a Type II Chebyshev IIR lowpass filter with the following specifications :
Passband frequency 1,200 Hz

Stopband frequency 1,700 Hz

Sampling frequency 6,000 Hz

Passband ripple 0.50 dB

Attach (i) Magnitude response (ii) Phase response (iii) Impulse Response (iv) Pole Zero plot

Task-3: Designing of an IIR Elliptic Filter

Pass band ripple 0.5 dB

Stop band attenuation 40 dB


Pass band frequency 800 Hz

Stop band frequency 1,000 Hz

Sampling frequency 4,000 Hz

a) Show magnitude and phase response in the same window.


b) Attach impulse response
c) Attach pole zero plot
d) Obtain filter information.
Task-4: IIR filter designing by Pole-Zero placement method

Design an IIR filter with the following specifications: Pass band at pi/2 and stop bands at 0 and pi.
Take = 8000 .

a) Attach pole zero plot

b) Attach magnitude response

c) Attach phase response

d) Attach impulse response

e) Attach filter coefficients


LAB SESSION 10(b)
OBJECTIVE:
Object of this lab is to design different IIR filter using FDA tool.

THEORY:
Filter Design and Analysis Tool (FDA Tool) is Graphic User Interface for designing and
analyzing filters. It is used to design FIR and IIR filters by entering the desired filter
specifications, or by importing filter from MATLAB workspace or by adding, moving or
deleting poles and zeros. After designing a filter, the response can be viewed and analyses in
other Graphic User Interface tool named Filter Visualization Tool (FV Tool) linked with FDA
Tool. The different types of responses that can be viewed are listed below:

Magnitude response
Phase response
Group delay
Phase delay
Impulse response
Step response
Pole-zero plot
Zero-phase plot

OPENING FDA TOOL WINDOW:


FDA Tool can be opened using command:
fdatool

Figure A
The different steps involved in designing a filter using FDA Tool can be listed as:

1. Selection of type of response required


2. Selection of type of filter design
3. Specifying the filter order
4. Entering the filter specifications
5. Entering the magnitude specifications

After providing the information listed above, filter can be designed and its response can be viewed and
analysed.

The complete description of the FDA Tool window and different steps required to design a filter
are elaborated below:

1. Selecting response type: The desired response type is selected from the list of available
options, i.e., lowpass, highpass, bandpass, bandstop, differentiation, multiband, peaking etc.

2. Type of design method: The design can be of FIR or IIR filter. Depending upon whether
FIR Or IIR filter design is selected, further options are available in the dropdown menu. In
IIR filter design, the different options available in dropdown menu are as given below:

Butterworth
Chebyshev type I
Chebyshev type II
Elliptic
Maximally flat
Least Pth-norm
Const least Pth-norm

Equirriple
Least square
Window
Const least squares
Complex equiripple
Least Pth norm
Constrained equiripple
Generalized equiripple
Constrained band equirriple
Interpolated FIR
The options available depend upon the selection of response type.

3. Filter order: Under this, two options available are as


User-defined order: Under this option user has to enter the order of the filter.
Minimum order: The other option available for selecting the filter order is
minimum order. It is calculated by system itself.

4. Filter specifications: Depending upon the response type and design method selected,
the graphical representation of generalized filter specifications appear in the display
region of FDA Tool. These specifications are ‘Frequency Specifications’ and
‘Magnitude Specification’.
These specifications are provided by the user, as per filter design requirement, in the
appropriate blocks.

5. Designing filter: After all the requisite information is entered, a filter can be designed by
clicking the ‘Design Filter’ button available at the bottom of the window. Filter |
coefficients are calculated and magnitude response appears in the display region.

(Note: ‘Design Filter’ button will be disabled once the filter coefficients are computed.
This button will be enabled again in case any changes are made in the filter specifications.)

6. Displaying filter responses: Once the filter coefficients are calcu

lated as per the specifications provided by the user, the display region will show
magnitude response of the designed filter. The other filter response characteristics
can be viewed in the display region or FV Tool. The response to be viewed can be
selected from the different icons displayed on the toolbar shown in Figure below.

Figure : Different Response Icons on the Toolbar

(NOTE: The different responses for display can also be selected from the
‘Analysis’ menu on menu bar.)

7. Current filter information: The information about the designed filter is given in
the ‘Current Filter Information’ region of FDA Tool window as shown
in Figure A The information provided is about the ‘structure’, ‘order’, ‘stability’ and
‘source’
Storing a filter
The designed filter is stored by clicking ‘Store Filter’ button in
the ‘Current Filter Information’ region.
Filter manager
The ‘Filter Manager’ button opens up a new Filter Manager window
(Figure B) showing the list of filters stored. This window also has
options as: Edit current filter, Cascade, Rename, Remove and FV Tool.

Figure B Filter Manager Window

To cascade two or more filters, highlight the designed filters and press ‘Cascade’ button.
A new cascaded filter is added to the ‘Filter Manager’.

8. Saving and opening filter design session:


The filter design session is saved as MAT-file and can be used later on. It can be saved by
clicking save icon or selecting save session
option in File menu and giving desired session name. Similarly, the saved session
can be opened by clicking open icon or by selecting open option in file menu and
selecting the previously saved filter design session.
FILTER VISUALIZATION TOOL:

The response characteristics can be viewed in a separate window by selecting the ‘Filter
Visualization Tool’ (FV Tool) from ‘view’ menu or clicking the ‘Full View Analysis’ button
on the toolbar. The FV Tool window is shown in Figure C

FV Tool has most of the menus on the menu bar and icons on the toolbar similar to that FDA
Tool with some additional icons which are mainly used to work with representation of the
responses.

Figure C Filter Visualization Tool (FV Tool) window


IIR FILTER DESIGN USING FDA TOOL

TASK-1
Design an IIR Butterworth band pass filter with the following specifications:

Normalized pass band edges at 0.40 and 0.65


Normalized stop band edges at 0.3 and 0.75
Pass band ripple 1 dB
Minimum stop band attenuation 40 dB

Show (i) Magnitude response (ii) Phase response (iii) Group delay (iv) Phase delay response.
Solution:
As per the given specifications, the requisite data is entered in new FDA Tool window
as shown in Figure

Figure. FDA Tool Window Showing Specification Entered and Magnitude Response for Task-1.

The filter is designed for minimum order so as to reduce the complexity of the design. In case, it
has to be designed for user defined order, then the order of the filter has to be calculated first by
user using appropriate formulas or MATLAB function.
The other responses can be viewed by clicking on the appropriate icon on the toolbar and
responses obtained are shown in Figures below

Figure Magnitude Response in dB


Figure Phase Response

Figure Group Delay Response


Figure. Phase Delay

These responses can also be viewed in FV Tool.

TASK-2
Design a Type II Chebyshev IIR lowpass filter with the following specifications :
Passband frequency 1,200 Hz

Stopband frequency 1,700 Hz

Sampling frequency 6,000 Hz


Passband ripple 0.50 dB

1. Show magnitude response.

2. Show pole/zero plot.

Solution: FDA Tool Window showing given specifications duly entered and magnitude
response in response display region is shown in Figure.

Figure. FDA Tool Window for Example


By using the FV Tool, the magnitude response and pole/zero plot is obtained as separate figures and is
shown in Figures.

Figure. Magnitude Response for Task-2.

Figure Pole/zero Plot for Q2


TASK-3:
Design an elliptic IIR low pass filter with following specifications :

Pass band ripple 0.5 dB

Stop band attenuation 40 dB

Pass band frequency 800 Hz

Stop band frequency 1,000 Hz

Sampling frequency 4,000 Hz


1. Show magnitude and phase response in the same window.
2. Obtain filter information.

Solution:
As per given specifications, the requisite data is entered in FDA Tool window. By
clicking appropriate icon on toolbar, the magnitude and phase responses are obtained in
the same window.
1. These magnitude and phase responses obtained are viewed in FV Tool window also
and are shown in Figure .
Figure . Magnitude and Phase Response for Task-3.

2. To obtain the information about the filter ‘Filter Information’ icon on Toolbar of FDA
Tool Window is clicked or ‘Filter Information’ option is selected from ‘Analysis’ menu. The
detailed filter information appears in the display region as shown in Figure a, b and c.
Figure ‘Filter information’ for Task-2
The filter information is obtained by scrolling down the text in the window
shown in Figure 15.41b.
EXERCISE:
Record Your Voice at home while turn any motor of your house ON.
Design a filter using FDA Tool.
Remove sound of motor from recorded signal.
Listen the output signal.

OUTCOME BASED EDUCATION


LAB RUBRICS
Digital Signal Processing (EE-394)
Student Name: Roll No: Section:

Excellent Good Average Poor Comments


100% 75% 50% 0%
The background The concepts The theoretical There is no
knowledge of
Clarity of theory is fully are not knowledge or
the background
Concepts grasped. completely concepts are
theory/concepts.
cleared. inappropriate.
The logic of
MATLAB the
program/Simulink program/model
The The
model is properly is correct but
Obtaining program/model program/model
simulated to there is some
plot/results has some is not correct.
obtain the desired syntax error to
major errors.
response. obtain the
results

Lab task Late


Submission submission of
submitted in time.
Lab Task

You might also like