DEPARTMENT OF ELECTRONICS AND COMMUNICATION

ENGINEERING
ANDHRA UNIVERSITY COLLEGE OF ENGINEERING
ANDHRA UNIVERSITY
VISAKHAPATNAM – 530003

DIGITAL SIGNAL PROCESSING LAB MANUAL

(B. Tech. 3rd Year – 2nd Semester)
(Effective from the admitted batch of 2021-2022)

1
 DIGITAL SIGNAL PROCESSING LAB
[EC 3207]
(Effective from the admitted batch of 2021-2022)

Periods Sessional Exam Total

Credits Exam Hrs.
Theory Tutorial Lab Marks Marks Marks

1.5 - - 3 3 50 50 100

Course Objectives:

• To make familiar with practical implementation of the digital signal processing.

• To develop DSP algorithms for convolution, correlation and DFT.
• To design digital filters.
• To have hands on experience in MATLAB and DSP processor.
Course Outcomes:

• Generation and Implementation of discrete time signals and systems using MATLAB
• Analyze the Frequency analysis of discrete signals and systems using MATLAB.
• Design and simulate FIR and IIR filters with different techniques using MATLAB.
• Verification of Linear and Circular Convolution using DSP Processor.
• Implementation of FIR and IIR filters with different techniques using DSP Processor.

SYLLABUS

1. Sampling theorem, illustration of up sampling in time and frequency domain

2. Sampling theorem, illustration of up sampling in time and frequency domain

3. Implement

a) Linear Convolution of two sequences

b) Circular Convolution of two sequences

c) Cross-Correlation and Auto-Correlation

4. FFT of a given N-point sequence using

0
 a) DIT-FFT b) DIF-FFT

5. System Response of discrete time sequences

a) Impulse b) Step

6. Spectral Analysis of given waveform and plot Spectrogram (Time V/s Frequency)

a) Sine b) Square c) Audio file

7. Study of architecture of DSP Chip TMS320C6711

8. Design following IIR Digital Filters using i) Butterworth and ii) Chebyshev designs:

a) LPF b) HPF c) BPF d) BSF

9. Design FIR digital filters using i) Rectangular Window ii) Hamming Window

a) LPF b) HPF c) BPF d) BSF

10. Addition of White Gaussian Noise to an audio file and recover the signal using
Butterworth filters.

11. Perform various operations on digital images

a) Cropping b) Rotation c) Histogram d) Binary image e) RGB to Gray Conversion

f) Water Marking g) Addition of noise to an image

1
 INDEX
S.No. Experiment Name of the Experiment Page No
No.
1 1 Sampling theorem 1-5
2 2 (a) Illustration of Up Sampling 6-7
3 2 (b) Illustration of Down Sampling 8-9
4 3 (a) Linear Convolution of two sequences 10-12
5 3 (b) Circular Convolution of two sequences 13-16
6 3 (c) Correlation of two sequences 17-19
7 4 (a) DIT-FFT 20-23
8 4 (b) DIF-FFT 24-27
9 5 (a) Impulse response of a system 28-30
10 5 (b) Step response of a system 31-33
11 6 (a) Spectrogram of Sine Wave 34-35
12 6 (b) Spectrogram of Square Wave 36-37
13 6 (c) Spectrogram of Audio file 38-39
14 7 Study of architecture of TMS320711 DSP chip 40-43
15 8 (a) IIR Butterworth Digital Filter 44-47
16 8 (b) IIR Chebyshev Digital Filter 48-52
17 9 FIR Digital Filter 53-58
18 10 Addition of Noise to an Audio File and recovery 59-62
of the Audio File using Butterworth Filter
19 11 (a) Image Cropping 63-64
20 11 (b) Image Rotation 65-66
21 11 (c) Histogram equalization 67-68
22 11 (d) Binary image 69-70
23 11 (e) RGB to Gray Conversion 71-73
24 11 (f) Addition of Noise to an image 74-75

2
3
4
 EXPERIMENT-1
SAMPLING THEOREM
Aim: To Verify the sampling theorem using MATLAB.

Software Required: MATLAB software 7.0 and above.

Theory:

Sampling Theorem: The sampling theorem specifies the minimum sampling rate at which a
continuous-time signal needs to be uniformly sampled so that the original signal can be
completely recovered or reconstructed by these samples alone.

fs ≥ 2fm

Nyquist rate: Nyquist rate states that the sampling signal frequency should be double the
input signal’s highest frequency component to get distortion less output signal. The theorem
states that for reconstructing a sampled signal accurately from the available samples, the
sampling frequency should be at least twice as much as the highest frequency component of
the signal.

fs = 2fm

Over sampling: It occurs if the sampling signal frequency is greater than double the input
signal’s highest frequency component.

fs > 2fm

Under sampling: The sampling signal frequency is less than double the input signal’s highest
frequency component, at this condition aliasing occurs. The spectra of X(ω) are overlapped in
this scenario because the sampling rate is less than the Nyquist rate, making it impossible to
extract the original signal from the sampled signal.

fs < 2fm

Applications of sampling:
• To maintain sound quality in music recordings.
• Sampling process applicable in the conversion of analog to discrete form.
• Speech recognition systems and pattern recognition systems.

1
 • Modulation and demodulation systems
• In sensor data evaluation systems
• Radar and radio navigation system sampling is applicable.

Procedure:
• Open MATLAB.
• Open New m-file by clicking on “New” and then “Script” in the Editor Window.
• Type the program and add comments to the program.
• Save the file with .m extension.
• Execute the program by clicking on “Run” in the editor window.
• Check the command window for errors and then enter the inputs as asked by the
command window.
• Observe the outputs in the command window and in Figure window.

MATLAB Code:
clc;
clear all;
close all;
t=-0.05:.000001:0.05; % input signal period
fm=input('Enter the frequency in Hz'); % taking the input in
Hz
x=cos(2*pi*fm*t);
subplot(4,4,1); %Taking an output plot with 4x4 figures
plot(t,x);
xlabel('time');
ylabel('x(t)');
title('continous time signal in Hz');
grid;
n1=-4:1:4;
fs1=1.6*fm;
fs2=2*fm;
fs3=8*fm;
x1=cos(2*pi*fm/fs1*n1); %Samples are spaced at a distance of
1/1.6 = 0.625
subplot(4,4,2);
stem(n1,x1);
xlabel('time');
ylabel('x(n)');
title('discrete time signal with fs<2fm');
hold on;
subplot(4,4,2);
plot(n1,x1);
grid;
n2=-5:1:5;
2
x2=cos(2*pi*fm/fs2*n2); % Samples are spaced at a distance of
1/2 = 0.5
subplot(4,4,3);
stem(n2,x2);
xlabel('time');
ylabel('x(n)');
title('discrete time signal with fs=2fm');
hold on;
subplot(4,4,3);
plot(n2,x2)
grid;
n3=-20:1:20;
x3=cos(2*pi*fm/fs3*n3); %Samples are taken at a distance of
1/8 = 0.125
subplot(4,4,4);
stem(n3,x3);
xlabel('time');
ylabel('x(n)');
title('discrete time signal with fs>2fm')
hold on;
subplot(4,4,4);
plot(n3,x3)
grid;
t2=-0.0005:.000001:0.0005; % Increase the time axis for -0.5ms
to 0.5ms
fm2=input('Enter the frequency in kHz-');
x4=cos(2*pi*fm2*1000*t2); % convert entered value into KHz by
multiplying with 1000
subplot(4,4,5);
plot(t2,x4);
xlabel('time');
ylabel('x(t)');
title('continous time signal in KHz');
grid;
n1=-4:1:4;
fs4=1.6*fm2;
fs5=2*fm2;
fs6=8*fm2;
x5=cos(2*pi*fm2/fs4*n1);
subplot(4,4,6);
stem(n1,x5);
xlabel('time');
ylabel('x(n)');
title('discrete time signal with fs<2fm');
hold on;
subplot(4,4,6);
plot(n1,x5);
grid;
n2=-5:1:5;
x6=cos(2*pi*fm2/fs5*n2);
subplot(4,4,7);
3
stem(n2,x6);
xlabel('time');
ylabel('x(n)');
title('discrete time signal with fs=2fm');
hold on;
subplot(4,4,7);
plot(n2,x6)
grid;
n3=-20:1:20;
x7=cos(2*pi*fm2/fs6*n3);
subplot(4,4,8);
stem(n3,x7);
xlabel('time');
ylabel('x(n)');
title('discrete time signal with fs>2fm');
hold on;
subplot(4,4,8);
plot(n3,x7)
grid;
t3=-0.0000005:.00000001:0.0000005; % Increase the time axis
for -0.5us to 0.5us
fm3=input('Enter the frequency in MHz-');
x8=cos(2*pi*fm3*1000000*t3); % convert entered value into KHz
by mulplying with 1000000
subplot(4,4,9);
plot(t3,x8);
xlabel('time');
ylabel('x(t)');
title('continous time signal in MHz');
grid;
n1=-4:1:4;
fs7=1.6*fm3;
fs8=2*fm3;
fs9=8*fm3;
x9=cos(2*pi*fm3/fs7*n1);
subplot(4,4,10);
stem(n1,x9);
xlabel('time');
ylabel('x(n)');
title('discrete time signal with fs<2fm');
hold on;
subplot(4,4,10);
plot(n1,x9);
grid;
n2=-5:1:5;
x10=cos(2*pi*fm3/fs8*n2);
subplot(4,4,11);
stem(n2,x10);
xlabel('time');
ylabel('x(n)');
title('discrete time signal with fs=2fm');
4
hold on;
subplot(4,4,11);
plot(n2,x10)
grid;
n3=-20:1:20;
x11=cos(2*pi*fm3/fs9*n3);
subplot(4,4,12);
stem(n3,x11);
xlabel('time');
ylabel('x(n)');
title('discrete time signal with fs>2fm')
hold on;
subplot(4,4,12);
plot(n3,x11);
grid;

OUTPUT:

i) Command Window:

Enter the frequency in Hz50

Enter the frequency in kHz-5

Enter the frequency in MHz-5

ii) Figure Window:

RESULT: Hence the Sampling Theorem was verified for 50 Hz, 5 KHz and 5 MHz signals
using MATLAB.

5
 EXPERIMENT – 2(a)
ILLUSTRATION OF UPSAMPLING
AIM:

To perform up-sampling on a given input signal using MATLAB

SOFTWARE REQUIRED:

MATLAB 7.0 and above

THEORY:

In many areas of digital signal processing (DSP) applications—such as communications,

speech, and audio processing—rising or lowering of a sampling rate is required. The
principle that deals with changing the sampling rate belongs essentially to mul rate signal
processing. The process of increasing sampling rate by an integer factor is known as
upsampling. This process is described as follows:
y(m) ={ x(m/L) m=NL
0 otherwise

where n = 0, 1, 2, … , x (n ) is the sequence to be upsampled by a factor of L , and y (m ) is

the upsampled sequence. As an example, suppose that the data sequence is given as follows:
x (n ):8 8 4 –5 –6 …
After upsampling the data sequence x (n ) by a factor of 3 (adding L – 1 zeros for each
sample), we have the upsampled data sequence w (m ) as:
w (m ): 8 0 0 8 0 0 4 0 0 –5 0 0 –6 0 0 …
The next step is to smooth the upsampled data sequence via an interpolation filter.

Procedure:
• Open MATLAB.
• Open New m-file by clicking on “New” and then “Script” in the Editor Window.
• Type the program and add comments to the program.
• Save the file with .m extension.
• Execute the program by clicking on “Run” in the editor window.
• Check the command window for errors and then enter the inputs as asked by the
command window.
• Observe the outputs in the command window and in Figure window.

MATLAB Code:
clc;
clear all;
close all;
6
ti
f = input('Enter the frequency of the input signal: '); %Take
the input signal frequency
sam = input('Enter the sample factor: '); %Take the input
sample factor
n= 0:(1/(f*10)):5*(1/f);
x=cos(2*pi*f*n);
y = zeros(1, sam*length(x)); %Generate 0's of sample fator
times time 'n'
y([1:sam:length(y)]) =x; %take x samples and ever factor
iteration
subplot(2,1,1);
stem(n,x);
title('input sequence');
xlabel('time index');
ylabel('amplitude');
subplot(2,1,2);
stem(n,y([1:length(x)]));
title('output sequence');
xlabel('time index');
ylabel('amplitude');

OUTPUTS:

Command window

Enter the frequency of the input signal: 1000

Enter the sample factor: 3

Figure Window:

input sequence
1

0.5
amplitude

-0.5

-1
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
time index 10 -3

output sequence
1

0.5
amplitude

-0.5

-1
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
time index 10 -3

RESULT: Hence the signal is upsampled successfully using MATLAB.

7
 EXPERIMENT – 2(b)
ILLUSTRATION OF DOWNSAMPLING
AIM:

To perform down-sampling on a given input signal using MATLAB

SOFTWARE REQUIRED:

MATLAB 7.0 and above

THEORY:

In many areas of digital signal processing (DSP) applications—such as communications,

speech, and audio processing—rising or lowering of a sampling rate is required. The
principle that deals with changing the sampling rate belongs essentially to multirate signal
processing.

The process of reducing a sampling rate by an integer factor is referred to

as downsampling of a data sequence. We also refer to downsampling as decimation (not
taking one of ten). The term decimation used for the downsampling process has been
accepted and used in many textbooks and fields. To downsample a data sequence x (n ) by an
integer factor of M , we use the following notation:
y (m ) = x (mM ),
where y (m ) is the downsampled sequence, obtained by taking a sample from the data
sequence x (n ) for every M samples (discarding M – 1 samples for every M samples). As an
example, if the original sequence with a sampling period T = 0.1 second (sampling rate = 10
samples per sec) is given by
x (n ):8 7 4 8 9 6 4 2 –2 –5 –7 –7 –6 –4 …
and we downsample the data sequence by a factor of 3, we obtain the downsampled sequence
as
y (m ):8 8 4 –5 –6 … ,
Procedure:
• Open MATLAB.
• Open New m-file by clicking on “New” and then “Script” in the Editor Window.
• Type the program and add comments to the program.
• Save the file with .m extension.
• Execute the program by clicking on “Run” in the editor window.
• Check the command window for errors and then enter the inputs as asked by the
command window.
8
 • Observe the outputs in the command window and in Figure window.

MATLAB Code:
clc;
close all;
clear all;
f = input('Enter the frequency of the input signal: ')
sam = input('Enter the sample factor: ')
n= 0:(1/(f*10)):5*(1/f);
x=cos(2*pi*f*n);
y=x([1:sam:length(x)]);
subplot(2,1,1);
stem(n,x);
title('input sequence');
xlabel('time index n');
ylabel('amplitude');
subplot(2,1,2);
stem(y);
title('output sequence');
xlabel('time index n');
ylabel('amplitude');

OUTPUTS:
a) Command window3
Enter the frequency of the input signal: 1000
Enter the sample factor: 3

b) Figure Window

input sequence
1

0.5
amplitude

-0.5

-1
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
time index n 10 -3

output sequence
1

0.5
amplitude

-0.5

-1
0 2 4 6 8 10 12 14 16 18
time index n

RESULT: Hence the signal is down-sampled successfully using MATLAB.

9
10
 EXPERIMENT – 3(a)
LINEAR CONVOLUTION OF TWO SEQUENCES
AIM:

To perform linear convolution of two sequences using MATLAB.

SOFTWARE REQUIRED:

MATLAB 7.0 and above

THEORY:

Linear convolution is a mathematical operation done to calculate the output of any Linear-
Time Invariant (LTI) system given its input and impulse response. We can represent Linear
Convolution as y(n)=x(n)*h(n). Here, y(n) is the output (also known as convolution sum).
x(n) is the input signal, and h(n) is the impulse response of the LTI system.

To calculate the number of samples in the output of linear convolution. [L = M + N – 1]

Where M is the number of samples in x(n). N is the number of samples in h(n).

Procedure:
• Open MATLAB.
• Open New m-file by clicking on “New” and then “Script” in the Editor Window.
• Type the program and add comments to the program.
• Save the file with .m extension.
• Execute the program by clicking on “Run” in the editor window.
• Check the command window for errors and then enter the inputs as asked by the
command window.
• Observe the outputs in the command window and in Figure window.

MATLAB Code:

clc;
clear all;
close all;
x1 = input('Enter the input sequence: ');
h1 = input('Enter the impulse response: ');
n1 = length(x1); %n1 is length of sequence x1
n2 = length(h1); %n2 is length of sequence h1
x2 = [x1, zeros(1,n1)]; % Add zeros to extend the sequence
length
h2 = [h1, zeros(1,n2)];
% To perform convolution
for i = 1:(n1+n2-1)
y(i) = 0;
11
for j = 1:n2
if (i-j+1)>0
y(i) = y(i) + x2(j)*h2(i-j+1);
end
end
end
disp('The convolution of x1 and h1 is:');
disp(y)
figure(1);
stem(x1);
title('Input sequence');
figure(2);
stem(h1);
title('Impulse response sequence');
figure(3);
stem(y);
title('output sequence');

OUTPUTS:

Command Window:

Enter the input sequence: [1 2 3 4]

Enter the impulse response: [1 2 1 2]

The convolution of x1 and h1 is:

1 4 8 14 15 10 8

Figure Window:

12
 Input sequence
4

3.5

2.5

1.5

0.5

0
1 1.5 2 2.5 3 3.5 4

Impulse response sequence

2

1.8

1.6

1.4

1.2

0.8

0.6

0.4

0.2

0
1 1.5 2 2.5 3 3.5 4

13
 output sequence
15

10

0
1 2 3 4 5 6 7

RESULT: Hence Linear Convolution is performed between two sequences successfully

using MATLAB.

14
 EXPERIMENT-3(b)
CIRCULAR CONVOLUTION OF TWO SEQUENCES
AIM:

To perform circular convolution of two sequences using MATLAB.

SOFTWARE REQUIRED:

MATLAB 7.0 and above

THEORY:

Circular convolution is essentially the same process as linear convolution. Just like linear
convolution, it involves folding a sequence, shifting it, multiplying it with another sequence,
and summing the resulting products. However, in circular convolution, the signals are all
periodic. Thus, the shifting can be thought of as actually being a rotation. Since the values
keep repeating because of the periodicity. Hence, it is known as circular convolution. Circular
Convolution is represented as y(n)=x(n)⊕h(n).

In circular convolution, both the sequences (input and impulse response) must be of equal
sizes. They must have the same number of samples. Thus the output of a circular convolution
has the same number of samples as the two inputs. If the sequences are not of equal length
they are made equal by inserting zeros in the sequence

Procedure:
• Open MATLAB.
• Open New m-file by clicking on “New” and then “Script” in the Editor Window.
• Type the program and add comments to the program.
• Save the file with .m extension.
• Execute the program by clicking on “Run” in the editor window.
• Check the command window for errors and then enter the inputs as asked by the
command window.
• Observe the outputs in the command window and in Figure window.

MATLAB Code:
clc;
clear all;
close all;
x = input('Enter the input sequence: ');
h = input('Enter the impulse response: ');
n1 = length(x); %n1 is length of sequence x
n2 = length(h); %n2 is length of sequence h
n = max(n1,n2); % maximum of n1 and n2 is considered as n

15
y = zeros(1,n); % Add zeros to extend the sequence length of
size n
if(n1<n2)
x = [x, zeros(1,n-n1)]; %if n1<n2 extend length of x
end
if(n1>n2)
h = [h, zeros(1, n-n2)]; %if n1>n2 extend length of y
disp(x);
disp(h);
end

for i=1:n
for j = 1 : n
z = mod(i-j, n);
y(i) = y(i) + x(j).*h(z+1);
disp(y);
end
end
figure(1);
stem(x);
title("Input sequence");
figure(2);
stem(h);
title("Impulse response sequence");
figure(3);
stem(y);
title("output sequence");

OUTPUTS:

Command Window:

Enter the input sequence: [1 2 3 4]

Enter the impulse response: [1 2 3]

1 2 3 4

1 2 3 0

18 16 10 16

16
Figure Window:

Input sequence
4

3.5

2.5

1.5

0.5

0
1 1.5 2 2.5 3 3.5 4

17
 Impulse response sequence
3

2.5

1.5

0.5

0
1 1.5 2 2.5 3 3.5 4

output sequence
18

16

14

12

10

0
1 1.5 2 2.5 3 3.5 4

RESULT: Hence the Circular Convolution is performed between two sequences using
MATLAB.
18
 EXPERIMENT-3(c)

CORRELATION OF TWO SEQUENCES

AIM:

To perform auto correlation and cross correlation between two sequences using MATLAB.

SOFTWARE REQUIRED:

MATLAB 7.0 and above

THEORY:

Correlation is a measure of the degree to which two sequences are similar.

Auto Correlation

Auto-correlation, sometimes known as serial correlation in the discrete time case, is the
correlation of a signal with a delayed copy of itself as a function of delay. Informally, it is the
similarity between observations of a random variable as a function of the time lag between
them. The analysis of autocorrelation is a mathematical tool for finding repeating patterns,
such as the presence of a periodic signal obscured by noise, or identifying the missing
fundamental frequency in a signal implied by its harmonic frequencies. It is often used in
signal processing for analyzing functions or series of values, such as time domain signals.
∞
r(x,x)(l) =
∑
x(n)x(n − l )
n=−∞

Cross Correlation

Cross-correlation mathematically measures the similarity of two different signals. Consider

an example where of having a set of data samples represented by x[n] and y[n]. Cross
correlation is used to measure on a sample-by-sample basis how similar x[n] is to y[n].
∞
r(x,x)(l) =
∑
x(n)y(n − l )
n=−∞

Procedure:
• Open MATLAB.
• Open New m-file by clicking on “New” and then “Script” in the Editor Window.
• Type the program and add comments to the program.
• Save the file with .m extension.
• Execute the program by clicking on “Run” in the editor window.
• Check the command window for errors and then enter the inputs as asked by the
command window.
19
 • Observe the outputs in the command window and in Figure window.

MATLAB Code:cclc;
close all;
clear all;
%Input sequence
x1=input('Enter x(n) sequence=');
x2=input('Enter h(n) sequence=');
N1=length(x1);
N2=length(x2);
n1=0:1:N1-1;
n2=0:1:N2-1;
y1=xcorr(x1,x2);%Cross Correlation
y2=xcorr(x1);%Auto Correlation
disp('Cross-Correlation of x(n) and h(n) sequences');
disp(y1);
disp('Auto-Correlation of x(n) sequence');
disp(y2);
%Plotting inputs and outputs
subplot(2,2,1);
stem(n1,x1);
xlabel('SEQUENCE');
ylabel('AMPLITUDE');
title('x(n) SEQUENCE');
subplot(2,2,2);
stem(n2,x2);
xlabel('SEQUENCE');
ylabel('AMPLITUDE');
title('h(n) SEQUENCE');
subplot(2,2,3);
stem(y1);
xlabel('SEQUENCE');
ylabel('AMPLITUDE');
title('CROSS CORRELATION OF x(n)& h(n)');
subplot(2,2,4);
stem(y2);
xlabel('SEQUENCE');
ylabel('AMPLITUDE');
title('AUTO CORRELATION OF x(n)');

OUTPUTS:

Command Window:

Enter x(n) sequence=[1 2 3 4]

Enter h(n) sequence=[1 2 1 2]

20
Cross-Correlation of x(n) and h(n) sequences

2.0000 5.0000 10.0000 16.0000 12.0000 11.0000 4.0000

Auto-Correlation of x(n) sequence

4.0000 11.0000 20.0000 30.0000 20.0000 11.0000 4.0000

Figure Window:

x(n) SEQUENCE h(n) SEQUENCE

4 2

3 1.5
AMPLITUDE

2 AMPLITUDE 1

1 0.5

0 0
0 1 2 3 0 1 2 3
SEQUENCE SEQUENCE
CROSS CORRELATION OF x(n)& h(n) AUTO CORRELATION OF x(n)
20 30

15
AMPLITUDE

AMPLITUDE

20
10
10
5

0 0
0 2 4 6 8 0 2 4 6 8
SEQUENCE SEQUENCE

RESULT: Hence auto correlation and cross correlation is performed between two sequences
successfully using MATLAB.

21
 EXPERIMENT – 4(a)
DIT-FFT (Decimation in Time – Fast Fourier Transform)
AIM:

To find the FFT of a given sequence through the DITFFT algorithm using MATLAB.

SOFTWARE REQUIRED:

MATLAB 7.0 and above

THEORY:

FFT: A fast Fourier transform (FFT) is an algorithm that computes the discrete Fourier
transform (DFT) of a given input sequence, or its inverse (IDFT). The complexity involved in
DFT is reduced through FFT via twiddle factor WNk considerations.

DFT is calculated using the equation

N−1
X(k) = x(n)e −j2πkn/N where k = 0 to N − 1
∑
n=0

The complex term is replaced with a value known as the twiddle factor. FFT equation is
given as follows
N−1
X(k) = x(n)WNk where k = 0 to N − 1
∑
n=0

FFT is calculated using two algorithms i) Decimation in Time (DIT) and ii) Decimation in
Frequency (DIF)

DIT-FFT:

• DITFFT algorithms are based upon the decomposition of the input sequence into smaller
and

smaller sub-sequences.

• In this input sequence x(n) is splitted into even and odd numbered samples.

• Splitting operation is done on time domain sequence.

• In DIT FFT input sequence is in bit reversed order while the output sequence is in natural

order.

22
The formula for calculating butterfly structures is as follows:

Upper-leg=A+B*WNk and

Lower-leg = A-B*WNk .

Procedure:
• Open MATLAB.
• Open New m-file by clicking on “New” and then “Script” in the Editor Window.
• Type the program and add comments to the program.
• Save the file with .m extension.
• Execute the program by clicking on “Run” in the editor window.
• Check the command window for errors and then enter the inputs as asked by the
command window.
• Observe the outputs in the command window and in Figure window.

MATLAB Code:

%****To find the DFT using DIT-FFT algorithm****%

clc; %clears command window
clear all; %clears the variables in the workspace
close all; %closes all the opened figure windows
tic; %tic-toc is used to calculate elapsed time for code
%execution. tic is mentioned at starting of the code and
toc %at the ending of the code%

% ****Input sequence****
xn=input('Enter the sequence: '); %Enter the input sequence for calculating DIT-FFT
N=input('Enter the length: '); %Enter the length of the sequence
p=length(xn); %Storing the length of the sequence in the variable p
if N>p %if the length of the sequence 'N' required is greater
than %the entered input sequence length 'p'
xn=[xn,zeros(1,N-p)]; %then append N-p zeros to 'xn'.
else
xn=xn; % otherwise 'xn' will remain same and return the
entered %sequence as it is
end %ending if loop
23
%*****Computing DFT using Radix-2 DITFFT (Decimation in Time FFT)
algorithm*****%
lev=log2(N); % No. of levels in DIT-FFT computation, Ex: for N=8,
%lev= log2(8) =3.
x=bitrevorder(xn); % input is considered in DITFFT in Bit-reversal order
tw=exp((-j*2*pi*(0:(N/2)-1))/N); %Twiddle factor used to calculate FFT
%**Calculating DIT-FFT algorithm using for loop**%
for levn=1:lev; %calculating FFT for 1:lev stages, Ex: for N=8,
calculating fft for 1 to 3 stages
L=2^levn; % defining the radix value 'L' for each stage of fft
calculation
twl=tw(1:N/L:N/2); %defining N/2 twiddle factors for N point FFT.
for k=0:L:N-L; %defining k range for calculating A and B values
for n=0:L/2-1; %defining n range for calculating A and B values
A=x(n+k+1); %defining upper leg input value 'A' in butterfly
structure
B=x(n+k+(L/2)+1)*twl(n+1); %defining lower leg input value 'B' in butterfly
structure
%Formulae for each stage of calculation in DIT-FFT using butterfly structure is
%upper level output = A+(twiddle factor)*B
%lower level output = A-(twiddle factor)*B
x(n+k+1)=A+B; %calculating upper level output of butterfly
x(n+k+(L/2)+1)=A-B; %calculating lower level output of butterfly
end
end
end %ending all the for loops used
x % Displaying the DIT-FFT value in command window
l=0:N-1; %Defining x-axis range for graph plotting

%****Plotting inputs and outputs of DIT-FFT****

%Four figures should be displayed. So, dividing figure window into 4 rows and 1 column
subplot(4,1,1); %plotting graph in 1st position of (4*1) figure window
stem(l,xn); %plotting input sequence 'xn'
title('Input Sequence'); %naming the graph with the title 'input sequence'
xlabel('n'); %labelling x axis as 'n'
ylabel('X[n]'); %labelling y axis as X[n]
subplot(4,1,2); %plotting graph at 2nd position of (4*1) figure window
stem(l,abs(x)); %plotting magnitude response of DIT-FFT i.e., of x
title('Magnitude Response'); %giving title of the graph as 'Magnitude response'
xlabel('k'); %labelling x axis as 'k'
ylabel('X[k]'); %labelling y axis as 'X[k]'
subplot(4,1,3); %plotting graph at 3rd position of (4*1) figure window
stem(angle(x)); %plotting phase response of DIT-FFT i.e., of x
title('Phase Response'); %naming the graph with the title 'Phase Response'
xlabel('k'); %labelling x axis as 'k'
24
ylabel('Phase(X[k])'); %labelling y axis as 'Phase X[k]'
subplot(4,1,4); %plotting graph at 4th position of (4*1) figure window
z=ifft(x,N); %calculating inverse FT of FFT sequence 'x' using ifft
disp('z=')
disp(z) %displaying inverse fft sequence in command window
stem(abs(z)); %plotting graph for inverse fft
title('Input by inverse fft'); %naming the graph with title 'input by inverse fft’
xlabel('n'); %labelling x axis as 'n'
ylabel('X[n]'); %labelling y axis as 'X[n]'
toc; %calculating time taken to execute the code

OUTPUTS:

Command Window:

Enter the sequence: [1 2 3 4]

Enter the length: 8

x=

Columns 1 through 5

10.0000 + 0.0000i -0.4142 - 7.2426i -2.0000 + 2.0000i 2.4142 - 1.2426i -2.0000 + 0.0000i

Columns 6 through 8

2.4142 + 1.2426i -2.0000 - 2.0000i -0.4142 + 7.2426i

z=

Columns 1 through 5

1.0000 + 0.0000i 2.0000 + 0.0000i 3.0000 + 0.0000i 4.0000 + 0.0000i 0.0000 + 0.0000i

Columns 6 through 8

0.0000 - 0.0000i 0.0000 - 0.0000i -0.0000 - 0.0000i

Figure window:

25
 Input Sequence
4

X[n]
2
0
0 1 2 3 4 5 6 7
n
Magnitude Response
10
X[k]
5
0
0 1 2 3 4 5 6 7
k
Phase Response
Phase(X[k])

2
0
-2
1 2 3 4 5 6 7 8
k
Input by inverse fft
4
X[n]

2
0
1 2 3 4 5 6 7 8
n

RESULT: Hence FFT for a given N-point sequence is computed through DIT-FFT algorithm
successfully using MATLAB.

26
 Experiment – 4(b)
DIF-FFT (Decimation in Frequency – Fast Fourier Transform)
AIM:

To find the FFT of a given sequence through the DIF-FFT algorithm using MATLAB.

SOFTWARE REQUIRED:

MATLAB 7.0 and above

THEORY:

FFT: A fast Fourier transform (FFT) is an algorithm that computes the discrete Fourier
transform (DFT) of a given input sequence, or its inverse (IDFT). The complexity involved in
DFT is reduced through FFT via twiddle factor WNk considerations.

DFT is calculated using the equation

N−1
X(k) = x(n)e −j2πkn/N where k = 0 to N − 1
∑
n=0

The complex term is replaced with a value known as the twiddle factor. FFT equation is
given as follows
N−1
X(k) = x(n)WNk where k = 0 to N − 1
∑
n=0

FFT is calculated using two algorithms i) Decimation in Time (DIT) and ii) Decimation in
Frequency (DIF)

DIF-FFT:
• DIFFFT algorithms are based upon the decomposition of the output sequence into
smaller andsmaller sub-sequences.
• In this output sequence X(k) is splitted into even and odd numbered samples
• Splitting operation is done on frequency domain sequence.
• In DIFFFT, input sequence is in natural order. And the output DFT should be
read in bit reversal order

27
The formula for calculating butterfly structures is as follows:

Upper-leg=A+B and

Lower-leg = (A-B)*WNk .

Procedure:
• Open MATLAB.
• Open New m-file by clicking on “New” and then “Script” in the Editor Window.
• Type the program and add comments to the program.
• Save the file with .m extension.
• Execute the program by clicking on “Run” in the editor window.
• Check the command window for errors and then enter the inputs as asked by the
command window.
• Observe the outputs in the command window and in Figure window.

MATLAB Code:

%To find the DFT using DIF-FFT

tic;
clc; %clears command window
clear all; %clears the variables in the workspace
close all; %closes all the opened figure windows

%Input sequence
X=input('Enter the sequence : '); %Enter the input sequence for calculating DIF-FFT
N=input('Enter the Point : '); %Enter the length of the sequence
n=length(X); %storing the length of the sequence X in the variable n
x=[X zeros(1,N-n)]; % appending zeros for satisfying lengths of N and X
M=log2(N); % No.of stages in DIF-FFT computation, Ex: for N=8,
%M=log2(8)=3.

%Computing DFT using DIF-FFT

for m=1:M %DIF-FFT is calculated for 1 to M stages, Ex: N=8
%sequence has 3 stages
d=2^(M-m+1); % determining radix-2 value for each stage,
%Ex: for N=8 sequence, i)stage 1,d=8 ii)stage 2,d=4
%iii)stage 3,d=2
for l=1:d:(N-d+1) %defining l value for DIF-FFT butterfly calculation
for k=0:(d/2)-1 %defining No.of twiddle factors for DIF-FFT
%k ranges from 0 to 3 for N=8
w=exp(-1i*2*pi*k/d); % Twiddle factor
z1=x(l+k); %determining upper-leg input for DIF-FFT butterfly
z2=x(l+k+(d/2)); %determining lower-leg input for DIF-FFT butterfly
%Formulae for each stage of calculation in DIF-FFT using butterfly structure is
28
%upper level output = z1+z2
%lower level output = (z1-z2)*twiddle factor
x(l+k)=z1+z2; %butterfly upper leg output
x(l+k+(d/2))=(z1-z2)*w; %butterfly lower leg output
end
end
end
y=bitrevorder(x); %for DIF-FFT, output is considered in bitreversal order
disp(y) % displaying the DIF-FFT output with ‘y’
%Plotting inputs and outputs
subplot(4,1,1) ; %plotting graph in 1st position of (4*1) figure window
stem(abs(X)); %plotting the input sequence
title('Input Sequence'); %naming the graph with title’Input Sequence’
xlabel('n'); %labelling X-axis as n
ylabel('X[n]'); %labelling Y-axis as ‘X[n]’
subplot(4,1,2); %plotting graph in 2nd position of (4*1) figure window
stem(abs(y)); %plotting the magnitude response of DIF-FFT
sequence
title('Magnitude Response');
xlabel('k'); %labelling X-axis as ‘k’
ylabel('X[k]'); %labelling Y-axis as ‘X[k]’
subplot(4,1,3); %plotting graph in 3rd position of (4*1) figure window
stem(angle(y)); %plotting the phase response of DIF-FFT sequence
title('Phase Response'); %naming the graph with title’Phase Response’
xlabel('k'); %labelling X-axis as ‘k’
ylabel('X[k]'); %labelling Y-axis as ‘X[k]’
subplot(4,1,4); %plotting graph in 4th position of (4*1) figure window
z=ifft(y,N) %Computing inverse fft
stem(abs(z)); %plotting the input sequence through ifft
title('Input by inverse fft'); %naming the graph with title’Input by inverse fft’
xlabel('n'); %labelling X-axis as ‘n’
ylabel('X[n]'); %labelling Y-axis as ‘X[n]’
toc; %calculating time taken to execute the code

OUTPUTS:
Command Window:
Enter the sequence : [1 2 3 4 1 2 3 4]
Enter the Point : 8
y=
Columns 1 through 5
20.0000 + 0.0000i 0.0000 + 0.0000i -4.0000 + 4.0000i 0.0000 + 0.0000i -4.0000 +
0.0000i
Columns 6 through 8 0.0000 + 0.0000i -4.0000 - 4.0000i 0.0000 + 0.0000i
z=
Columns 1 through 5
1.0000 + 0.0000i 2.0000 - 0.0000i 3.0000 + 0.0000i 4.0000 + 0.0000i 1.0000 + 0.0000i
Columns 6 through 8
29
2.0000 - 0.0000i 3.0000 + 0.0000i 4.0000 + 0.0000i

Figure Window:

Input Sequence
4
X[n]

2
0
1 2 3 4 5 6 7 8
n
Magnitude Response
20
X[k]

10
0
1 2 3 4 5 6 7 8
k
Phase Response
2
X[k]

0
-2
1 2 3 4 5 6 7 8
k
Input by inverse fft
4
X[n]

2
0
1 2 3 4 5 6 7 8
n

RESULT: Hence FFT for a given N-point sequence is calculated via DIF-FFT algorithm
successfully using MATLAB.

30
 EXPERIMENT - 5 (a)
IMPULSE RESPONSE OF A SYSTEM
AIM :

To find the impulse response of a system using MATLAB.

SOFTWARE REQUIRED :

MATLAB 7.0 and above

THEORY:

A discrete time system performs an operation on an input signal based on a predefined

criteria to produce a modified output signal. The input signal x(n) is the system excitation and
y(n) is the system response.

If the input to the system is a unit impulse i.e., x(n)=δ(n) then the output of the system is
known as impulse response denoted by h(n) where h(n)=T[δ(n)].

Procedure:
• Open MATLAB.
• Open New m-file by clicking on “New” and then “Script” in the Editor Window.
• Type the program and add comments to the program.
• Save the file with .m extension.
• Execute the program by clicking on “Run” in the editor window.
• Check the command window for errors and then enter the inputs as asked by the
command window.
• Observe the outputs in the command window and in Figure window.

MATLAB Code:
clc;
close all;
clear all;
a=input('Enter the coefficients of x(n),x(n-1),.... =');
%enter the coefficients of x(n)
b=input('Enter the coefficients of y(n),y(n-1),.... =');
%enter the coefficients of y(n)
N=input('Enter the Number of Samples=');%Enter number of
samples required
t=(0:N-1)*1 %Time Vector
syms z h(n); %creating symbolic variables z and h(n)
H = poly2sym(b,(z))/poly2sym(a,(z))%Transfer function of H(z)
h(n) = iztrans(H) %Inverse Z-Transform
h=h(t) % evaluating h at each point of t
c=double(h) %converting values of h into double precision
31
%plotting inputs and outputs
%plotting
subplot(3,1,1); %plotting h(n) at 1st position in 3*1 figure
window
stem(t,h); %plotting discrete sequence h
title('h(n)'); % naming the graph as h(n)
xlabel('n'); %labelling x-axis as n
ylabel('h(n)'); %labelling y-axis as h(n)
%Unit_Impulse
subplot(3,1,2);%plotting unit impulse at 2nd position in 3*1
figure window
x=[ones(1,1),zeros(1,N-1)] %defining unit impulse function
stem(t,x); %plotting unit-impulse in discrete form
title('Unit Impulse'); % naming the graph as unit impulse
xlabel('n'); %labelling x-axis as n
ylabel('x(n)'); %labelling y-axis as x(n)
%Convoluting h(n) and x(n) to obtain y(n)
subplot(3,1,3); %plotting y(n) at 3rd position in 3*1 figure
window
y=conv(c,x) %convolution of h(n) and x(n)
t1=(0:length(y)-1)*1 %defining time vector of y(n)
stem(t1,y); %plotting discrete y(n)
title('y(n)=h(n)*x(n)'); %naming the graph
xlabel('n'); %labelling x-axis as n
ylabel('y(n)'); %labelling y-axis as y(n)

OUTPUTS:

Command Window:
Enter the coefficients of x(n),x(n-1),.... =[1 2 3 4]

Enter the coefficients of y(n),y(n-1),.... =[1 2 3]

Enter the Number of Samples=10

t=

0 1 2 3 4 5 6 7 8 9

H=

(z^2 + 2*z + 3)/(z^3 + 2*z^2 + 3*z + 4)

h(n) =

(3*kroneckerDelta(n, 0))/4 + symsum(-(2*r3*r3^n + 3*r3^n + r3^n*r3^2)/(2*r3^2 + 6*r3 +

12), r3 in root(z1^3 + 2*z1^2 + 3*z1 + 4, z1))

h=

[ 0, 1, 0, 0, -4, 8, -4, 0, -20, 56]

c=
32
0 1 0 0 -4 8 -4 0 -20 56

x=

1 0 0 0 0 0 0 0 0 0

y=

Columns 1 through 17

0 1 0 0 -4 8 -4 0 -20 56 0 0 0 0 0 0 0

Columns 18 through 19

0 0

t1 =

Columns 1 through 17

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

Columns 18 through 19

17 18

Figure Window:

h(n)
40
h(n)

20
0
-20
0 1 2 3 4 5 6 7 8 9
n
Unit Impulse
1
x(n)

0.5

0
0 1 2 3 4 5 6 7 8 9
n
y(n)=h(n)*x(n)
40
y(n)

20
0
-20
0 2 4 6 8 10 12 14 16 18
n

33
RESULT: Hence the impulse response of a system is evaluated successfully using MATLAB.

34
 EXPERIMENT - 5 (b)
STEP RESPONSE OF A SYSTEM
AIM:

To find the step response of a system using MATLAB.

SOFTWARE REQUIRED:

MATLAB 7.0 and above

THEORY:

A discrete time system performs an operation on an input signal based on a predefined

criteria to produce a modified output signal. The input signal x(n) is the system excitation and
y(n) is the system response.

If the input to the system is a unit step i.e., x(n)=u(n) then the output of the system is known
as step response. It represents how a system behaves when subjected to a sudden change in
input from zero to one. Mathematically, the step response of a system is obtained by
convolving the system's impulse response with a unit step function. The unit step function,
denoted as u(n), is defined as u(n)=0 if n<0 and 1 if n≥0.

Procedure:
• Open MATLAB.
• Open New m-file by clicking on “New” and then “Script” in the Editor Window.
• Type the program and add comments to the program.
• Save the file with .m extension.
• Execute the program by clicking on “Run” in the editor window.
• Check the command window for errors and then enter the inputs as asked by the
command window.
• Observe the outputs in the command window and in Figure window.

MATLAB Code:

clc;
close all;
clear all;
%Accepting inputs from user
a=input('Enter the coefficients of x(n),x(n-1),.... =');
b=input('Enter the coefficients of y(n),y(n-1),.... =');
N=input('Enter the Number of Samples=');
t=(0:N-1)*1%Time Vector
syms z h(n);
H = poly2sym(b,(z))/poly2sym(a,(z))%Transfer function of H(z)
h(n) = iztrans(H)%Inverse Z-Transform
35
h=h(t)
c=double(h)
%h(n)
subplot(3,1,1);
stem(t,h);
title('h(n)');
xlabel('n');
ylabel('h(n)');
%Unit_Step
subplot(3,1,2);
x=[ones(1,N)]
stem(t,x);
title('Unit Step');
xlabel('n');
ylabel('x(n)');
%Convolution
subplot(3,1,3);
y=conv(c,x)
t1=(0:length(y)-1)*1
stem(t1,y);
title('y(n)=h(n)*x(n)');
xlabel('n');
ylabel('y(n)');

OUTPUTS:

Command Window:

Enter the coefficients of x(n),x(n-1),.... =[1 3 5 7]

Enter the coefficients of y(n),y(n-1),.... =[2 4 6 8]
Enter the Number of Samples=10
t=
0 1 2 3 4 5 6 7 8 9
H=
(2*z^3 + 4*z^2 + 6*z + 8)/(z^3 + 3*z^2 + 5*z + 7)
h(n) =
(8*kroneckerDelta(n, 0))/7 - symsum(-(4*r3*r3^n + 6*r3^n + 2*r3^n*r3^2)/(3*r3^2 + 10*r3
+ 21), r3 in root(z1^3 + 3*z1^2 + 5*z1 + 7, z1))
h=
[ 2, -2, 2, -2, 10, -34, 66, -98, 202, -578]
c=
2 -2 2 -2 10 -34 66 -98 202 -578

36
x=
1 1 1 1 1 1 1 1 1 1
y=
Columns 1 through 17
2 0 2 0 10 -24 42 -56 146 -432 -434 -432 -434 -432 -442 -408 -474
Columns 18 through 19
-376 -578
t1 =
Columns 1 through 17
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
Columns 18 through 19
17 18
Figure Window:

h(n)
200
0
h(n)

-200
-400

0 1 2 3 4 5 6 7 8 9
n
Unit Step
1
x(n)

0.5

0
0 1 2 3 4 5 6 7 8 9
n
y(n)=h(n)*x(n)
0
y(n)

-200
-400

0 2 4 6 8 10 12 14 16 18
n

RESULT: Hence the step response of a system is evaluated successfully using MATLAB.

37
 EXPERIMENT – 6(a)
SPECTROGRAM OF SINE WAVE
AIM:

To plot spectrogram of a sine wave using MATLAB.

SOFTWARE USED:

MATLAB 7.0 and above

THEORY:

A spectrogram is a visual way of representing the signal strength, or “loudness”, of a signal

over time at various frequencies present in a particular waveform. Not only can one see
whether there is more or less energy at, but one can also see how energy levels vary over
time.

In other sciences spectrograms are commonly used to display frequencies of sound waves
produced by humans, machinery, animals, whales, jets, etc., as recorded by microphones.

In the seismic world, spectrograms are increasingly being used to look at frequency content
of continuous signals recorded by individual or groups of seismometers to help distinguish
and characterize different types of earthquakes or other vibrations in the earth.

Procedure:
• Open MATLAB.
• Open New m-file by clicking on “New” and then “Script” in the Editor Window.
• Type the program and add comments to the program.
• Save the file with .m extension.
• Execute the program by clicking on “Run” in the editor window.
• Check the command window for errors and then enter the inputs as asked by the
command window.
• Observe the outputs in the command window and in Figure window.

MATLAB Code:
clc;
clear all;
close all;
%Accepting the input values from the user
fs= input('Enter the sampling frequency=');
f= input('Enter the frequency of sine waveform=');
N=input('Enter the number of samples=');
t =(0:(N-1))/fs;%Time vector
x=sin(2*pi*f*t);%Sinusoid Signal
38
spectrogram(x,[],[],[],fs,'yaxis');%Plotting the spectrogram
colorbar;%Displaying the colorbar

OUTPUTS:

Command Window:

Enter the sampling frequency=10000

Enter the frequency of sine waveform=2500

Enter the number of samples=2000

Figure Window:

RESULT: Hence the spectrogram of Sine wave is plotted using MATLAB.

39
 EXPERIMENT – 6(b)
SPECTROGRAM OF SQUARE WAVE
AIM:

To plot spectrogram of a square wave using MATLAB.

SOFTWARE USED:

MATLAB 7.0 and above

THEORY:

A spectrogram is a visual way of representing the signal strength, or “loudness”, of a signal

over time at various frequencies present in a particular waveform. Not only can one see
whether there is more or less energy at, but one can also see how energy levels vary over
time.

In other sciences spectrograms are commonly used to display frequencies of sound waves
produced by humans, machinery, animals, whales, jets, etc., as recorded by microphones.

In the seismic world, spectrograms are increasingly being used to look at frequency content
of continuous signals recorded by individual or groups of seismometers to help distinguish
and characterize different types of earthquakes or other vibrations in the earth.

Procedure:
• Open MATLAB.
• Open New m-file by clicking on “New” and then “Script” in the Editor Window.
• Type the program and add comments to the program.
• Save the file with .m extension.
• Execute the program by clicking on “Run” in the editor window.
• Check the command window for errors and then enter the inputs as asked by the
command window.
• Observe the outputs in the command window and in Figure window.

MATLAB Code:
clc;
clear all;
close all;
%Accepting the input values from the user
fs= input('Enter the sampling frequency=');
f= input('Enter the frequency of sine waveform=');
N=input('Enter the number of samples=');
t =(0:(N-1))/fs;%Time vector
x=square(2*pi*f*t);%Square Signal
40
spectrogram(x,[],[],[],fs,'yaxis');%Plotting the spectrogram
colorbar;%Displaying the colorbar

OUTPUTS:

Command Window:

Enter the sampling frequency=10000

Enter the frequency of sine waveform=2500

Enter the number of samples=2000

Figure Window:

RESULT: Hence the spectrogram of square wave is plotted using MATLAB.

41
 EXPERIMENT – 6(c)
SPECTROGRAM OF AUDIO FILE
AIM:

To plot the spectrogram of an audio file using MATLAB.

SOFTWARE USED:

MATLAB 7.0 and above

THEORY:

A spectrogram is a visual way of representing the signal strength, or “loudness”, of a signal

over time at various frequencies present in a particular waveform. Not only can one see
whether there is more or less energy at, but one can also see how energy levels vary over
time.

In other sciences spectrograms are commonly used to display frequencies of sound waves
produced by humans, machinery, animals, whales, jets, etc., as recorded by microphones.

In the seismic world, spectrograms are increasingly being used to look at frequency content
of continuous signals recorded by individual or groups of seismometers to help distinguish
and characterize different types of earthquakes or other vibrations in the earth.

Procedure:
• Open MATLAB.
• Open New m-file by clicking on “New” and then “Script” in the Editor Window.
• Type the program and add comments to the program.
• Save the file with .m extension.
• Execute the program by clicking on “Run” in the editor window.
• Check the command window for errors and then enter the inputs as asked by the
command window.
• Observe the outputs in the command window and in Figure window.

MATLAB Code:
clc;
clear all;
close all;
[x,fs] = audioread('male.wav');%Reading the wav file from
specified location
spectrogram(x,[],[],[],fs,'yaxis');%Plotting the spectrogram
colorbar;%Displaying the colorbar

42
OUTPUT:

4
-40

3.5 -50

-60
3
-70
Frequency (kHz)

2.5
-80

2 -90

-100
1.5
-110
1
-120

0.5 -130

-140
0
5 10 15 20 25 30 35 40 45
Time (secs)
RESULT: Hence the spectrogram of an audio file is plotted using MATLAB.

43
 EXPERIMENT -7
STUDY OF THE ARCHITECTURE OF TMS320C6711 DSP CHIP

AIM: To study the architecture of DSP TMS320C6711 chip.

Differences between DSP and Microprocessor:

• A microprocessor incorporates the functions of a computer's central processing unit
(CPU) on a single or few integrated circuits (IC). The purpose of a microprocessor is
to accept digital data as input, process it as per the instructions, and then provide the
output. The DSP processor, on the other hand, is a particular type of microprocessor.
It processes a digital signal or information signal.
• The main difference between a DSP and a microprocessor is that a DSP has features
designed to support high performance, repetitive, numerically intensive tasks.
• Microprocessor is often used for computation, text editing, multimedia display, and
communication over a network while DSP is used in applications such as image
processing, speech recognition and telecommunications.
• As compared to general microprocessors, DSP processors are more efficient at
performing basic arithmetic operations, especially multiplication

Features of DSP:
• Parallel multiply and add.
• Multiple memory access (to fetch two operands and store the result).
• Lots of registers to hold data temporarily.
• DSP aims to modify or improve the signal. It is characterized by the representation of
discrete units, such as discrete time, discrete frequency, or discrete domain signals.
• Ability to move data to and from memory quickly.

Introduction:
• The Texas Instruments TMS320C6000 generation of high-performance DSPs now
includes the TMS320C6711.
• The ‘C6711 floating point device will begin sampling in fourth quarter 1999,
providing 900 MFLOPS (million floating-point operations per second) at 150MHz,
with a lower cost version available at 100MHz.
• This processor provides 8 execution units, including 2 multipliers and 6 arithmetic
logic units (ALUs).
• These units operate in parallel and can perform up to 6 floating-point instructions
during a single clock cycle up to 900 MFLOPS at 150MHz
44
• This is a 32-bit floating point DSP.
• On chip memory-72k bytes.

COMPONENTS

EXTERNAL MEMORY INTERFACE:

• This 32-bit interface supports direct interface with external devices like SBSRAM,
SDRAM, SRAM, ROM.
• Supports interface up to 4 banks of 64M bit SDRAM.
• Also supports 8-bit, 16-bit and 32-bit external interfaces in asynchronous mode.

HOST PORT INTERFACE:

• Provides direct access to memory mapped peripherals of the TMS320C6711
processor.
• One DMA channel is dedicated to it.

SERIAL PORTS:
• This processor has 2 serial ports.
• Fully independent transmit and receive sections.
• Time Division Multiplexed (TDM) mode.
45
TIMERS:
• Two 32-bit timers used to time and count events or to interrupt the CPU
• directs an external ADC to start conversion or the DMA controller to start a data
transfer.
• includes a time period register, which specifies the timer’s frequency; a timer counter
register, which contains the value of the incrementing counter; and a timer control
register, which monitors the timer’s status.

PIPELINING:

Used to get parallel instructions working properly, requiring careful timing

There are three stages of pipelining: program fetch, decode, and execute.

1. The program fetch stage is composed of four phases:

a) PG: program address generate (in the CPU) to fetch an address

b) PS: program address send (to memory) to send the address.

c) PW: program address ready wait (memory read) to wait for data.

d) PR: program fetch packet receive (at the CPU) to read opcode from memory.

2.The decode phase is composed of two phases

a) DP: to dispatch all the instructions within an FP to the appropriate functional units.

b) DC: instruction decode.

3. The execute stage is composed of 6 phases (with fixed point) to 10 phases (with floating
point), due to delays (latencies) associated with the following instructions:

a) Multiply instruction, which consists of 2 phases due to 1 delay.

b) Load instruction, which consists of 5 phases due to 4 delays.

c) Branch instruction, which consists of 6 phases due to 5 delays.

REGISTER FILES:
• The computational core consists of two register files.
• Provides local storage for arithmetic input data, results.
• Each register file consists of sixteen 32 bit registers.

FUNCTIONAL UNITS:
• CPU consists of 8 units divided into 2 data paths A, B.
• (.M) for multiply operations.
• (.L) for logical and arithmetic operations.
46
 • (.S) for branch, bit manipulation, and arithmetic operations.
• (.D) for loading/storing and arithmetic operations.

Each functional unit can read directly from or write directly to the register file within its own
path.
• Each path includes a set of sixteen 32-bit registers, A0 through A15 and B0 through
B15.
• Units ending in 1 write to register file A, and units ending in 2 write to register file B.
• Each functional unit side can access data from the registers on the opposite side using
a cross-path.

TYPES OF INSTRUCTIONS:
• Add/Subtract/Multiply

• Branch/move

47
RESULT: Hence the architecture of TM320C6711 DSP Chip is studied.ss

48
 EXPERIMENT – 8(a)
IIR DIGITAL BUTTERWORTH FILTER
AIM: To design IIR digital Butterworth LPF, HPF, BPF, and BSF filters using MATLAB.

SOFTWARE REQUIRED: MATLAB 7.0 and above

THEORY:

A filter rejects unwanted frequencies from the input signals and allows the desired
frequencies. The filters are of different types
➢ Low-pass filter: Passes signals with frequencies lower than the cutoff frequency

➢ High-pass filter: Passes signals with frequencies higher than the cutoff frequency

➢ Band-pass filter: Passes signals within a specified frequency range

➢ Band-stop filter: Blocks signals within a specified range

A digital filter is a linear time-invariant discrete system.

Infinite Impulse Response (IIR) filters are of recursive type, whereby the present output
sample depends on the present input, past input samples, and output samples.

The different types of IIR filters are i) Butterworth filter ii) Chebyshev filter

Butterworth Filter:

The signal processing filter which has a flat frequency response in the passband can be
termed a Butterworth filter and is also called a maximally flat magnitude filter. The order of
the filter is given as

100.1αs − 1
log
100.1αp − 1
N ≥ Ω
log Ω s
p

Where

αs =stop band attenuation

αp= pass band attenuation
Ωs= pass band frequency in rad/sec
Ωp= stop band frequency in rad/sec
MATLAB Commands for Butterworth filter design:

[n,Wn] = buttord(Wp,Ws,Rp,Rs) returns the lowest order, n, of the digital Butterworth filter
with no more than Rp dB of passband ripple and at least Rs dB of attenuation in the
stopband. Wp and Ws are respectively the passband and stopband edge frequencies of the
49
filter, normalized from 0 to 1, where 1 corresponds to π rad/sample. The scalar (or vector) of
corresponding cutoff frequencies, Wn, is also returned. To design a Butterworth filter, use the
output arguments n and Wn as inputs to butter.

[b,a] = butter(n,Wn) returns the transfer function coefficients of an nth-order lowpass digital
Butterworth filter with normalized cutoff frequency Wn.

[b,a] = butter(n,Wn,ftype) designs a lowpass, highpass, bandpass, or bandstop Butterworth

filter, depending on the value of ftype and the number of elements of Wn.

Procedure:
• Open MATLAB.
• Open New m-file by clicking on “New” and then “Script” in the Editor Window.
• Type the program and add comments to the program.
• Save the file with .m extension.
• Execute the program by clicking on “Run” in the editor window.
• Check the command window for errors and then enter the inputs as asked by the
command window.
• Observe the outputs in the command window and in Figure window.

MATLAB Code:
clc;
clear all;
close all;
%Taking Input parameters
Ap=input('Enter passband attenuation in db=');
As=input('Enter stopband attenuation in db=');
fp=input('Enter passband frequency in hz=');
fs=input('Enter stopband frequency in hz=');
Fs=input('Enter sampling frequency in hz=');
w1=2*fp/Fs; %Passband corner frequency Wp
w2=2*fs/Fs; %Stopband corner frequency Ws
%Designing Butterworth filter
%LPF
[N,wn]=buttord(w1,w2,Ap,As) % N returns the lowest order
and wn returns the cutoff frequency of Butterworth filter
[b,a]=butter(N,wn); % [b,a] returns the transfer function
coefficients of an nth-order lowpass Butterworth filter
[h,o]=freqz(b,a,256); % [h,o] returns the frequency response
h and corresponding frequencies o
m=20*log10(abs(h)); %Convert magnitude into dB
subplot(2,2,1);
plot((o*Fs)/pi,m); % plotting frequency response
grid;
title('LPF');
50
ylabel('Gain in dB-->');
xlabel('(a)frequency-->');
%HPF
[N,wn]=buttord(w1,w2,Ap,As) % N returns the lowest order and
wn returns the cutoff frequency of Butterworth filter
[b,a]=butter(N,wn,'high'); % [b,a] returns the transfer
function coefficients of an nth-order lowpass Butterworth
filter
[h,o]=freqz(b,a,256); % [h,o] returns the frequency response
h and corresponding frequencies o
m=20*log10(abs(h)); %Convert magnitude into dB
subplot(2,2,2);
plot((o*Fs)/pi,m); % plotting frequency
response
grid;
title('HPF');
ylabel('Gain in dB-->');
xlabel('(b)frequency-->');
%BPF
[N,wn]=buttord(w1,w2,Ap,As); % N returns the lowest order
and wn returns the cutoff frequency of Butterworth filter
[b,a]=butter(N,0.5*[w1 w2]); % [b,a] returns the transfer
function coefficients of an nth-order lowpass Butterworth
filter
[h,o]=freqz(b,a,256); % [h,o] returns the frequency response
h and corresponding frequencies o
m=20*log10(abs(h)); %Convert magnitude into dB
subplot(2,2,3);
plot((o*Fs)/pi,m); % plotting frequency response
grid;
title('BPF');
ylabel('Gain in dB-->');
xlabel('(c)frequency-->');
%BSF
[N,wn]=buttord(w1,w2,Ap,As); % N returns the lowest order
and wn returns the cutoff frequency of Butterworth filter
[b,a]=butter(N,0.5*[w1 w2], 'stop'); % [b,a] returns the
transfer function coefficients of an nth-order lowpass
Butterworth filter
[h,o]=freqz(b,a,256); % [h,o] returns the frequency response
h and corresponding frequencies o
m=20*log10(abs(h)); %Convert magnitude into dB
subplot(2,2,4); % plotting frequency response
plot((o*Fs)/pi,m);
grid;
title('BSF');
ylabel('Gain in dB-->');
xlabel('(d)frequency-->');

OUTPUTS:

51
Command Window:

Enter passband attenuation in db=3

Enter stopband attenuation in db=30

Enter passband frequency in hz=1000

Enter stopband frequency in hz=2000

Enter sampling frequency in hz=5000

=3

wn

= 0.4914

Figure Window:

LPF HPF
0 0
Gain in dB-->

Gain in dB-->

-50 -50

-100 -100

-150 -150
0 2000 4000 6000 0 2000 4000 6000
(a)frequency--> (b)frequency-->
BPF BSF
0 0
Gain in dB-->

Gain in dB-->

-50
-50

-100

-100
-150
0 2000 4000 6000 0 2000 4000 6000
(c)frequency--> (d)frequency-->

RESULT: Hence IIR digital Butterworth LPF, HPF, BPF, BSF filters are designed using
MATLAB.

52
 EXPERIMENT – 8(b)

IIR DIGITAL CHEBYSHEV FILTER

AIM: To design IIR digital Chebyshev LPF, HPF, BPF, and BSF filters using MATLAB.

SOFTWARE REQUIRED: MATLAB 7.0 and above

THEORY:

A filter rejects unwanted frequencies from the input signals and allows the desired
frequencies. The filters are of different types
➢ Low-pass filter: Passes signals with frequencies lower than the cutoff frequency

➢ High-pass filter: Passes signals with frequencies higher than the cutoff frequency

➢ Band-pass filter: Passes signals within a specified frequency range

➢ Band-stop filter: Blocks signals within a specified range

A digital filter is a linear time-invariant discrete system.

Infinite Impulse Response (IIR) filters are of recursive type, whereby the present output
sample depends on the present input, past input samples, and output samples.

The different types of IIR filters are i) Butterworth filter ii) Chebyshev filter

Chebyshev filter:

There are two types of Chebyshev filters.Type 1 Chebyshev filters are all pole filters that
exhibit equiripple behavior in the passband and monotonic characteristics in the stopband.
Type II Chebyshev filters contain both poles and zeros and exhibit monotonic behavior in the
pass band and equiripple behavior in the stop band.

The order of the filter is given as

100.1αs − 1
cosh −1
100.1αp − 1
N ≥ Ω
cosh −1 Ω s
p

Where

αs =stop band attenuation

αp= pass band attenuation
Ωs= pass band frequency in rad/sec
Ωp= stop band frequency in rad/sec
MATLAB commands for Chebyshev filter

53
[n,Wp] = cheb1ord(Wp,Ws,Rp,Rs) returns the lowest order n of the Chebyshev Type I filter
that loses no more than Rp dB in the passband and has at least Rs dB of attenuation in the
stopband. The scalar (or vector) of corresponding cutoff frequencies Wp, is also returned.
Use the output arguments n and Wp with the cheby1 function.

[b,a] = cheby1(n,Rp,Wp) returns the transfer function coefficients of an nth-order lowpass

digital Chebyshev Type I filter with normalized passband edge frequency Wp and Rp
decibels of peak-to-peak passband ripple.

[b,a] = cheby1(n,Rp,Wp,ftype) designs a lowpass, highpass, bandpass, or bandstop

Chebyshev Type I filter, depending on the value of ftype and the number of elements of Wp.

Procedure:
• Open MATLAB.
• Open New m-file by clicking on “New” and then “Script” in the Editor Window.
• Type the program and add comments to the program.
• Save the file with .m extension.
• Execute the program by clicking on “Run” in the editor window.
• Check the command window for errors and then enter the inputs as asked by the
command window.
• Observe the outputs in the command window and in Figure window.

MATLAB Code:

clc;

clear all;

close all;

%Taking Input parameters

Ap=input('Enter passband attenuation in db=');

As=input('Enter stopband attenuation in db=');

fp=input('Enter passband frequency in hz=');

fs=input('Enter stopband frequency in hz=');

Fs=input('Enter sampling frequency in hz=');

w1=2*fp/Fs; %Passband corner frequency Wp

w2=2*fs/Fs; %Stopband corner frequency Ws

%Designing Chebyshev filter

%LPF
54
 [N,wn]=cheb1ord(w1,w2,Ap,As) % N returns the lowest order and wn returns the cutoff
frequency of Chebyshev Type I filter

[b,a]=cheby1(N,Ap,wn); %[b,a] returns the transfer function coefficients of an

nth-order lowpass Chebyshev Type I filter

[h,o]=freqz(b,a,256); % [h,o] returns the frequency response h and

corresponding frequencies o

m=20*log10(abs(h)); %Convert magnitude into dB

subplot(2,2,1);

plot((o*Fs)/pi,m); %plotting frequency response

grid;

title('LPF');

ylabel('Gain in dB-->');

xlabel('(a)frequency-->');

%HPF

[N,wn]=cheb1ord(w1,w2,Ap,As) % N returns the lowest order and wn returns the cutoff

frequency of Chebyshev Type I filter

[b,a]=cheby1(N,Ap,wn,’high’); %[b,a] returns the transfer function coefficients of an

nth-order lowpass Chebyshev Type I filter

[h,o]=freqz(b,a,256); % [h,o] returns the frequency response h and

corresponding frequencies o

m=20*log10(abs(h)); %Convert magnitude into dB

subplot(2,2,2);

plot((o*Fs)/pi,m); %plotting frequency response

grid;

title('HPF');

ylabel('Gain in dB-->');

xlabel('(b)frequency-->');

%BPF

[N,wn]=cheb1ord(w1,w2,Ap,As) % N returns the lowest order and wn returns the cutoff

frequency of Chebyshev Type I filter

55
[b,a]=cheby1(N,Ap,0.5*[w1 w2]); %[b,a] returns the transfer function coefficients of an
nth-order lowpass Chebyshev Type I filter

[h,o]=freqz(b,a,256); % [h,o] returns the frequency response h and

corresponding frequencies o

m=20*log10(abs(h)); %Convert magnitude into dB

subplot(2,2,3);

plot((o*Fs)/pi,m); %plotting frequency response

grid;

title('BPF');

ylabel('Gain in dB-->');

xlabel('(c) frequency-->');

%BSF

[N,wn]=cheb1ord(w1,w2,Ap,As) % N returns the lowest order and wn returns the cutoff

frequency of Chebyshev Type I filter

[b,a]=cheby1(N,Ap,0.5*[w1 w2],’stop’); %[b,a] returns the transfer function coefficients

of an nth-order lowpass Chebyshev Type I filter

[h,o]=freqz(b,a,256); % [h,o] returns the frequency response h and

corresponding frequencies o

m=20*log10(abs(h)); %Convert magnitude into dB

subplot(2,2,4);

plot((o*Fs)/pi,m); %plotting frequency response

grid;

title('BSF');

ylabel('Gain in dB-->');

xlabel('(d) frequency-->');

OUTPUTS:

Command Window:

Enter passband attenuation in db=3

Enter stopband attenuation in db=30

Enter passband frequency in hz=1000

56
Enter stopband frequency in hz=2000

Enter sampling frequency in hz=10000

N=

wn =

0.2000

Figure Window:

LPF HPF
0 0

-50
Gain in dB-->

Gain in dB-->

-50
-100

-150
-100
-200
0 5000 10000 0 5000 10000
(a)frequency--> (b)frequency-->
BPF BSF
0 0

-50
Gain in dB-->

Gain in dB-->

-100 -50

-150

-100
-200
0 5000 10000 0 5000 10000
(c) frequency--> (d) frequency-->

RESULT: Hence IIR digital Chebyshev LPF, HPF, BPF, BSF filters are designed using
MATLAB.

57
 EXPERIMENT – 9
FIR DIGITAL FILTERS
AIM: To design FIR digital LPF, HPF, BPF, and BSF filters using Rectangular and Hamming
Window using MATLAB.

SOFTWARE REQUIRED: MATLAB 7.0 and above

THEORY:

The FIR filters are of non-recursive type where by present output sample depends on present
input and previous input samples.

Based on frequency response, the filters can be classified as

➢ Low-pass filter: Passes signals with frequencies lower than the cutoff frequency

➢ High-pass filter: Passes signals with frequencies higher than the cutoff frequency

➢ Band-pass filter: Passes signals within a specified frequency range

➢ Band-stop filter: Blocks signals within a specified range

The advantages of FIR filter are:

➢ FIR filters have exact linear phase

➢ FIR filters are always stable

➢ FIR filters can be realized in both recursive and non-recursive structure

The three well-known methods for designing FIR filters with linear phase are i) Window
method ii) Frequency Sampling method iii) optimal or maximal design.

One of the most frequently used methods is the Window method. The different types of
windows are Rectangular, Kaiser, Hamming, and Hanning, Triangular, and Bartlett
Windowing techniques.

b = fir1(n,Wn,ftype) designs a lowpass, highpass, bandpass, bandstop, or multiband filter,

depending on the value of ftype and the number of elements of Wn.

PROCEDURE:
➢ Open MATLAB

➢ Click on the “New” in editor panel → select “Script”.

➢ Now Type the Program and also add comments to the program.

➢ Now save the file with a file name along with “.m” extension.

➢ Now Click on “Run” button present in the Editor panel.

➢ Give the inputs in the command window.

58
 ➢ Now observe the obtained waveforms in the figure window.

MATLAB PROGRAM:

%fir filt design window techniques

clc;

clear all;

close all;

rp=input('enter passband ripple in db=');

rs=input('enter the stopband ripple in db=');

fp=input('enter passband freq in hz=');

fs=input('enter stopband freq in hz=');

f=input('enter sampling freq in hz=');

wp=2*fp/f; %Passband corner frequency Wp

ws=2*fs/f; %Stopband corner frequency Ws

num=-20*log10(sqrt(rp*rs))-13;

dem=14.6*(fs-fp)/f;

n=ceil(num/dem);

n1=n+1;

if(rem(n,2)~=0)

n1=n;

n=n-1;

end

c=input('enter your choice of window function 1. rectangular 2. hamming : \n ');

if(c==1)

y=rectwin(n1);

disp('Rectangular window filter response');

end

if (c==2)

y=hamming(n1);

disp('hamming window filter response');

end
59
%LPF

b=fir1(n,wp,y); % b returns the transfer function coefficients of

an nth-order lowpass fir filter

[h,o]=freqz(b,1,256); % [h,o] returns the frequency response h and

corresponding frequencies o

m=20*log10(abs(h)); %converting magnitude into db

subplot(2,2,1);

plot((o*f)/pi,m);

grid;

title('LPF');

ylabel('Gain in dB-->');

xlabel('(a)frequency-->');

%HPF

b=fir1(n,wp,'high',y); % b returns the transfer function coefficients of

an nth-order highpass fir filter

[h,o]=freqz(b,1,256); % [h,o] returns the frequency response h and

corresponding frequencies o

m=20*log10(abs(h)); %converting magnitude into db

subplot(2,2,2);

plot((o*f)/pi,m);

grid;

title('HPF');

ylabel('Gain in dB-->');

xlabel('(b)frequency-->');

%BPF

wn=(1/2)*[wp ws];

b=fir1(n,wn,y); % b returns the transfer function coefficients of

an nth-order bandpass fir filter

[h,o]=freqz(b,1,256); % [h,o] returns the frequency response h and

corresponding frequencies o

m=20*log10(abs(h)); %converting magnitude into db

subplot(2,2,3);
60
plot((o*f)/pi,m);

grid;

title('BPF');

ylabel('Gain in dB-->');

xlabel('(c)frequency-->');

%BSF

b=fir1(n,wn,'stop',y); % b returns the transfer function coefficients of

an nth-order bandstop fir filter

[h,o]=freqz(b,1,256); % [h,o] returns the frequency response h and

corresponding frequencies o

m=20*log10(abs(h)); %converting magnitude into db

subplot(2,2,4);

plot((o*f)/pi,m);

grid;

title('BSF');

ylabel('Gain in dB-->');

xlabel('(d)frequency-->');

OUTPUTS:

1. RECTANGULAR WINDOW

Command Window:

enter passband ripple in db=0.1

enter the stopband ripple in db=0.001

enter passband freq in hz=3000

enter stopband freq in hz=6000

enter sampling freq in hz=10000

enter your choice of window function 1. rectangular 2. hamming :

Rectangular window filter response

61
Figure Window:

LPF HPF
0 0
Gain in dB-->

Gain in dB-->

-20 -20

-40 -40

-60
-60
0 5000 10000 0 5000 10000
(a)frequency--> (b)frequency-->
BPF BSF
0 0

-20 -5
Gain in dB-->

Gain in dB-->

-10
-40
-15
-60
-20
0 5000 10000 0 5000 10000
(c)frequency--> (d)frequency-->

2. HAMMING WINDOW

Command Window:

enter passband ripple in db=0.01

enter the stopband ripple in db=0.001

enter passband freq in hz=2000

enter stopband freq in hz=3000

enter sampling freq in hz=8000

62
enter your choice of window function 1. rectangular 2. hamming :

hamming window filter response

Figure Window:

LPF HPF
0 0
Gain in dB-->

Gain in dB-->

-50 -50

-100 -100
0 2000 4000 6000 8000 0 2000 4000 6000 8000
(a)frequency--> (b)frequency-->
BPF BSF
0 0
Gain in dB-->

Gain in dB-->

-50 -5

-100 -10
0 2000 4000 6000 8000 0 2000 4000 6000 8000
(c)frequency--> (d)frequency-->

RESULT: Hence FIR digital LPHF, HPF, BPF, and BSF filters are designed using
Rectangular and Hamming Window techniques in MATLAB.

63
 EXPERIMENT – 10
ADDITION OF NOISE TO AN AUDIO FILE AND RECOVERY OF THE
AUDIO FILE USING BUTTERWORTH FILTER
AIM: To add noise to an audio file and recover the signal using Butterworth filter using
MATLAB.

SOFTWARE REQUIRED: MATLAB 7.0 and above

THEORY:

This experiment involves adding white Gaussian noise to an audio file and recovering the
audio file using a digital bandpass Butterworth filter.

The steps involved are:

➢ Reading the audio file “male.wav” from the system

➢ Estimating power spectral density of audio file

➢ Addition of white Gaussian noise to the audio file

➢ Creating an audio file of the noisy signal

➢ Estimating power spectral density of noisy signal

➢ Implementing Butterworth bandpass filter to remove noise

➢ Determining lower and higher cutoff frequencies

➢ Filtering the noisy signal using a designed IIR filter

➢ Creating an audio file of the filtered audio file

➢ Estimating power spectral density of the filtered audio file.

[y,fs] = audioread(filename) reads data from the file named filename, and returns sampled
data, y, and a sample rate for that data, Fs. In this program, for “male.wav”, fs=8000 and
number of smaples is 408226.

PROCEDURE:
➢ Open MATLAB

➢ Click on the “New” in editor panel → select “Script”.

➢ Now Type the Program and also add comments to the program.

➢ Now save the file with a file name along with “.m” extension.

➢ Now Click on “Run” button present in the Editor panel.

➢ Give the inputs in the command window.

➢ Now observe the obtained waveforms in the figure window.

64
MATLAB PROGRAM:

clc;

close all;

clear all;

%Audio read

[y1,fs] = audioread('D:\DSP LAB for B.Tech \8\male.wav'); %Reading the wav file from
specified location

N=length(y1);

t =(0:(N-1))/fs; %Time vector

Y = fft(y1,N);

Pyy = Y.*conj(Y)/N; %Estimation of PSD of audio file

M=ceil(N/2);

f = fs/N*(0:M); %Frequency vector

subplot(6,1,1);

plot(f,Pyy(1:M+1));

axis([0 4000 0 1]);

title('Spectral plot of the signal')

xlabel('Frequency (Hz)');

ylabel('Power spectral density(db/Hz)');

subplot(6,1,2);

plot(t,y1);

title('Time Plot of the signal')

xlabel('Time(sec)');

ylabel('Amplitude');

%Addition of noise

x=awgn(y1,20); %Addition of white gaussian noise with an SNR of 20db

%Creating an audio file of the noisy signal

audiowrite ('D:\DSP LAB for B.Tech \8 \noise1.wav',x,fs);

X = fft(x,N); % FFT of noisy signal

Pxx = X.*conj(X)/N; %Estimation of PSD of signal+noise

65
subplot(6,1,3);

plot(f,Pxx(1:M+1));

axis([0 4000 0 1]);

title('Spectral plot of the signal+noise')

xlabel('Frequency (Hz)');

ylabel('Power spectral density(db/Hz)');

subplot(6,1,4);

plot(t,x);

title('Time Plot of the signal+noise')

xlabel('Time(sec)');

ylabel('Amplitude');

%Filter Design

n = 6; %Order of the filter

beginFreq = 300 / (fs/2); %Lower cut-off frequency of bandpass filter

endFreq =1100 / (fs/2); %Higher cut-off frequency of bandpass filter

%Calculating the digital bandpass butterworth filter

[b,a] = butter(n, [beginFreq,endFreq],'bandpass');

fOut = filter(b, a, x); %Filtering the noisy signal using the designed filter

%Creating an audio file of the filtered signal

audiowrite('D:\DSP LAB for B.Tech\8\result2.wav',fOut,fs);

F = fft(fOut,N);

Pff = F.*conj(F)/N; %Estimation of PSD of filtered signal

subplot(6,1,5);

plot(f,Pff(1:M+1));

axis([0 4000 0 1]);

title('Spectral plot of the filtered signal');

xlabel('Frequency (Hz)');

ylabel('Power spectral density(db/Hz)');

subplot(6,1,6);

plot(t,fOut);
66
title('Time Plot of the filtered signal')

xlabel('Time(sec)');

ylabel('Amplitude');

OUTPUT:

RESULT: Hence White gaussian noise is added to the audio file and the audio file is
recovered successfully using Bandpass butterworth filter using MATLAB.

67
 Experiment 11(a)
IMAGE CROPPING
Aim: To crop an image using MATLAB.
Software required: MATLAB 7.0 and above – Image Processing Toolbox
Procedure:

1. Start the MATLAB program.

2. Open new M-file

3. Type the program

4. Save in current directory

5. Compile and Run the program

6. If any error occurs in the program correct the error and run it again

7. For the output see command window\ Figure window

8. Stop the program.

MATLAB Program:
clc;
clear
all;close
all;
I=imread('circuit.tif'); %reading the image
I1=imcrop(I); %cropping the image using imcrop function
subplot(1,2,1)
imshow(I) %displaying original image
title('Original
Image')subplot(1,2,2)

68
imshow(I1) %dsiplaying cropped image
title('Cropped Image')
% Cropping image by specifying the crop rectangle as a four-element position vector,
%[xmin ymin width height]
I2=imcrop(I, [60 40 100 90]);
figure
imshow(I2) %dsiplaying cropped image
title('Cropped Image')

Output:

Result: Hence the image is cropped using MATLAB

69
 Experiment - 11 (b)
IMAGE ROTATION
Aim: To rotate an image at different angles using MATLAB.
Software required: MATLAB 7.0 and above – Image Processing Toolbox
Procedure:
1. Open MATLAB.

2. Open new M-file.

3. Type the program.

4. Save in current directory.

5. Compile and Run the program.

6. If any error occurs in the program correct the error and run it again.

7. For the output see command window\ Figure window.

8. Stop the program.

MATLAB Program:

clc;
clear
all;close
all;
img=imread('llama.jpg'); %reading the image
subplot(2,2,1);
imshow(img); %displaying original image
title('Original image');
I1=imrotate(img,45); %rotating the image at 45 degrees counterclockwise
subplot(2,2,2);
imshow(I1); %displaying the 45 degree rotated image
title('45 degree rotated image');
I2=imrotate(img,-90); %rotating the image at 90 degrees clockwise

70
subplot(2,2,3);
imshow(I2); %displaying the 90 degree rotated image
title('90 degree rotated image');
I3=imrotate(img,180); %rotating the image at 180 degrees counterclockwise
subplot(2,2,4);
imshow(I3); %displaying the 180-degree rotated image
title('180 degree rotated image');

Output:

Result: Hence the image is rotated at 45, 90, and 180 degrees using MATLAB.

71
 Experiment - 11(c)
HISTOGRAM EQUALIZATION
Aim: To perform histogram equalization using MATLAB.

Software required: MATLAB 7.0 and above – Image Processing Toolbox

Procedure:
1. Open MATLAB.

2. Open new M-file.

3. Type the program.

4. Save in current directory.

5. Compile and Run the program.

6. If any error occurs in the program correct the error and run it again.

7. For the output see command window\ Figure window.

8. Stop the program.

MATLAB

Program:clc;

clear

all;close

all;

I=imread('pout.tif'); %reads the image

subplot(2,2,1);

imshow(I); %displays original image

title('Original

image')

subplot(2,2,2);
imhist(I,64); %displays histogram of original image
title('histogram of original image')

J=histeq(I); %histogram equalization

subplot(2,2,3);

72
title('equalized image')

imshow(J); % displays histogram equalized image

subplot(2,2,4);

imhist(J,64); %displays histogram of equalized image

title('histogram of equalized image')

Output:

Result: Hence, histogram equalization is performed for a given image using MATLAB.

73
 Experiment -
11(d)BINARY
IMAGE
Aim: To convert an image into a binary image using MATLAB.
Software required: MATLAB 7.0 and above – Image Processing Toolbox

MATLAB Program:
%Binarize image using global threshold

I=imread('coins.png'); %reads the image

BW=imbinarize(I); %converting original image into binary image

figure

imshowpair(I,BW,'montage'); %displaying original and binary images side by side

title('Original & Binary images');

74
%Binarize Images with Darker Foreground than

BackgroundI = imread('printedtext.png'); %reads the

image figure
imshow(I) %displays original image
title('Original Image')

%converting into binary image

BW = imbinarize(I,'adaptive','ForegroundPolarity','dark','Sensitivity',0.4);

figure
imshow(BW) %diplays original image
title('Binary Version of Image')

Outputs:

Result: Hence an image is converted into a binary image using MATLAB.

75
 Experiment - 11(e)
RGB to GRAY CONVERSION
Aim: To convert RGB image into gray image using MATLAB.

Procedure:
1. Open MATLAB.
2. Open new M-file.
3. Type the program.
4. Save in current directory.
5. Compile and Run the program.
6. If any error occurs in the program correct the error and run it again.
7. For the output see command window\ Figure window.
8. Stop the program.

76
MATLAB Program:
clc; clear
all;close
all;
i=imread('onion.png'); %reads the image
subplot(3,2,1);
imshow(i); %displays the original RGB image
title('original RGB image');
%To display red channel, green and blue channels must be zero
a=i;
i(:,:,2)=0; %equating green channel to zero
i(:,:,3)=0; %equating blue channel to zero
subplot(3,2,2);
imshow(i); %display red channel
title('Red Channel');
%To display green channel, red and blue channels must be zero
i=a;
i(:,:,1)=0; %equating red channel to zero
i(:,:,3)=0; %equating blue channel to zero
subplot(3,2,3);
imshow(i); %display green channel
title('Green Channel');
%To display blue channel, red and green channels must be zero
i=a;
i(:,:,1)=0; %equating red channel to zero
i(:,:,2)=0; %equating green channel to zero
subplot(3,2,4);
imshow(i); %display blue channel
title('Blue Channel');

77
%RGB to gray conversion
j=imread('onion.png'); %read the image
figure
subplot(1,2,1)
imshow(j) % display original image
title('RGB image');
y=rgb2gray(j) %RGB to Gray conversion
subplot(1,2,2)
imshow(y) %displays gray image
title('Gray image')
Output:

Result: Hence Red, Green, and Blue channels in an image are displayed and RGB image is
converted into a Gray image.

78
 Experiment - 11(f)
ADDITION OF NOISE TO AN IMAGE
Aim: To add noise to the image and to recover the noisy image using mean and median
filtering using MATLAB.

Software required: MATLAB 7.0 and above – Image Processing Toolbox

1. MATLAB

Program:clc;
clear all;
close all;

a=imread('onion.png'); %reads the image

i1=rgb2gray(a); %RGB to GRAY conversion

subplot(2,2,1);

imshow(i1); % displays gary image

title('original image')

i2=imnoise(i1,'salt & pepper', 0.20); %addition of salt and pepper noise to an image

subplot(2,2,2);

imshow(i2); %displays noisy image

title('Noisy Image')

b=(1/9)*[1 1 1;1 1 1;1 1 1]; %mask for mean filter

i3=imfilter(i2,b); %mean filtering on noise image

subplot(2,2,3);

imshow(i3); %dsiplays mean filtered image

title('Filtering with mean filter')

i4=medfilt2(i2); %median filtering on noise image

subplot(2,2,4);

imshow(i4); %displays median fltered image

title('Filtering with median filter')

Output:

79
80