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

DSP Lab Manual final

Digital Signals

Uploaded by

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

DSP Lab Manual final

Digital Signals

Uploaded by

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

DEPARTMENT OF ELECTRONICS AND

COMMUNICATION ENGINEERING

LAB MANUAL

1
LIST OF EXPERIMENTS
CYCLE - 1 ............................................................................................................................ 2
GENERATION OF DISCRETE SIGNALS .....................................................................................................3
LINEAR CONVOLUTION ........................................................................................................................... 14
CIRCULAR CONVOLUTION ...................................................................................................................... 21
DISCRETE FOURIER TRANSFORM AND INVERSE DISCRETE FOURIER TRANSFORM................. 29
ADDITION OF SINUSOIDAL SIGNALS ................................................................................................... 36

CYCLE - 2 .................................................................................................................................................. 43
DESIGN OF ANALOG BUTTERWORTH FILTER .................................................................................... 44
DESIGN OF ANALOG CHEBYSHEV FILTER ............................................................................................ 52
DESIGN OF IIR DIGITAL FILTERS ............................................................................................................ 64
DESIGN OF FIR DIGITAL FILTERS USING RECTANGLE WINDOW ................................................... 82
DESIGN OF FIR DIGITAL FILTERS USING TRIANGLE WINDOW ....................................................... 90

CYCLE - 3……………………………………………………………………………………………………………………………………………………99
HISTOGRAM EQUALIZATION …….. ..................................................................................................... 100
EDGE DETECTION………………………………………………………………………………………………………………………..107

2
1
CYCLE - 1

2
GENERATION OF DISCRETE SIGNALS
AIM: To generate a discrete time signal (Impulse, Unit Step, Unit Ramp, Exponential, Parabola, Sinusoidal,
Square, Rectangle, Triangular and Sawtooth) using MATLAB.

SOFTWARE REQUIRED: MATLAB

THEORY:
Real signal can be quite complicated. The study of signals therefore starts with analysis of basic and
fundamental signals. For linear system, a complicated signal and its behaviour can be studied by superposition of
basic signals. Common basic signals are as follows

Discrete - Time signals:

Unit impulse sequence:


1; 𝐹𝑜𝑟 𝑛 = 0
𝛿(𝑛) = {
0; 𝑂𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
Unit step sequence:
1; 𝐹𝑜𝑟 𝑛 ≥ 0
𝑢(𝑛) = {
0; 𝑂𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
Unit ramp sequence:
𝑛; 𝐹𝑜𝑟 𝑛 ≥ 0
𝑟(𝑛) = {
0; 𝑂𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒

Exponential sequence:
𝑥 (𝑛 ) = 𝑎 𝑛

Sinusoidal sequence:
𝑥(𝑛) = 𝐴 sin(𝜔𝑛 + ф)
Square sequence:
𝑁
1, 0≤𝑛≤
𝑥 (𝑛 ) = { 2
𝑁
−1, <𝑛≤𝑁
2

Continuous - Time signals:

Unit impulse signals:


1; 𝐹𝑜𝑟 𝑡 = 0
𝛿(𝑡) = {
0; 𝑂𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
Unit step signals:
1; 𝐹𝑜𝑟 𝑡 ≥ 0
𝑢(𝑡) = {
0; 𝑂𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
Unit ramp signals:
𝑡; 𝐹𝑜𝑟 𝑡 ≥ 0
𝑟(𝑡) = {
0; 𝑂𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒
Exponential signals:
𝑥 (𝑡) = 𝑒 𝑎𝑡

3
Sinusoidal signals:
𝑥(𝑡) = 𝐴 sin(𝜔𝑡 + ф)
Square sequence:
𝑇
1, 0≤𝑡<
𝑥 (𝑡 ) = { 2
𝑇
−1, <𝑡≤𝑇
2
ALORITHM:
Step1: Start the program.
Step2: Get the input for signal generation.
Step3: Generate the waveform by using appropriate library function or logic.
Step4: Display the signal.

PROCEDURE:
1. Start the MATLAB program.
2. Create a NEW M-File from the file Menu.
3. An editor window open, type the program in it.
4. Save the file 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.

SOURCE CODE:
%Generation of discrete-time sequences
Clc;
clear all;
close all;
%% Impulse Function
n = input('enter value of n for Impulse Function: ');
x = -n:1:n;
y = [zeros(1,n),1,zeros(1,n)];
figure;
subplot(1,2,1)
plot(x,y),xlabel('t'), ylabel('x(t)')
grid on
title('Impulse signal');
subplot(1,2,2)
stem(x,y),xlabel('n'), ylabel('x(n)')
grid on
title('Impulse sequence');
%% Step Function
n = input('enter value of n for step function: ');
x = -n:1:n;

4
y = [zeros(1,n) ones(1,n+1)];
figure;
subplot (1,2,1);
plot (x,y),xlabel('t'), ylabel('u(t)'),axis([-1 n 0 2])
grid on
title('Unit step signal');
subplot (1,2,2);
stem (x,y),xlabel('n'), ylabel('u(n)'),axis([-1 n 0 2])
grid on
title('Unit step sequence');
%% Ramp Function
n = input('enter value of n for ramp function: ');
x = -n:1:n;
y = [zeros(1,n) 0:n];
figure;
subplot (1,2,1);
plot (x,y),xlabel('t'), ylabel('r(t)'),axis([-1 n 0 n])
title('Unit Ramp signal');
grid on
subplot (1,2,2);
stem (x,y),xlabel('n'), ylabel('r(n)'),axis([-1 n 0 n])
grid on
title('Unit Ramp sequence');
%% Exponential function
n = input('enter value of n for exponential function: ');
x = 0:0.5:n;
y = exp(x);
figure;
subplot (1,2,1);
plot(x,y),xlabel('t'), ylabel('x(t)')
title('Exponential signal');
grid on
subplot (1,2,2);
stem(x,y),xlabel('n'), ylabel('x(n)')
title('Exponential sequence');
grid on
%% Parabola function
n = input('enter value of n for parabola function: ');
x = -n:1:n;
y = x.^2;
figure;
subplot (1,2,1);
plot (x,y),xlabel('t'), ylabel('r(t)'),axis([-n n 0 (n^2)+1])
title('Parabolic signal');
grid on
subplot (1,2,2);

5
stem (x,y),xlabel('n'), ylabel('r(n)'),axis([-n n 0 (n^2)+1])
title('Parabolic sequence');
grid on
%% Sinusoidal function
x = 0:pi/100:2*pi;
y = sin(x);
figure;
subplot (1,2,1);
plot (x,y),xlabel('t'), ylabel('r(t)')
title('Sinusoidal signal');
grid on
subplot (1,2,2);
stem (x,y),xlabel('n'), ylabel('r(n)')
title('Sinusoidal sequence');
grid on
%% Square Function
n = input('enter time of plot for square function: ');
x = 1:1:n;
y = square(x);
figure
subplot (2,1,1);
plot (x,y),xlabel('t'), ylabel('x(t)')
title('Square signal');
grid on
subplot (2,1,2);
stem (x,y),xlabel('n'), ylabel('x(n)')
title('Square sequence');
grid on
%% Rectangular Function
duty = input('Enter duty cycle for rectangle function: ');
A = input('Enter amplitude for rectangle function: ');
a = 2*(duty/25);
b = 2*((100-duty)/25);
x = A*[-1*ones(1,b),ones(1,a)];
y = [x,x,x];
figure
subplot (2,1,1);
plot (y),xlabel('t'), ylabel('x(t)')
title('Rectangular signal');
grid on
subplot (2,1,2);
stem (y),xlabel('n'), ylabel('x(n)')
title('Rectangular sequence');
grid on
%% Triangular function
n = input('Enter amplitude of the Triangular waveform: ');

6
b = [0:n-1,n:-1:1];
a = [b,b,b,b];
figure
subplot (1,2,1);
plot (a),xlabel('t'), ylabel('x(t)')
title('Triangular signal');
grid on
subplot (1,2,2);
stem (a),xlabel('n'), ylabel('x(n)')
title('Triangular sequence');
grid on
%% Sawtooth Function
clc
clear all
close all
n = input('Enter amplitude of the sawtooth waveform: ');
x = 0:1:n;
y = x;
a = [y,y,y,y];
figure;
subplot (1,2,1);
plot (a),xlabel('t'), ylabel('r(t)')
title('Sawtooth signal');
grid on
subplot (1,2,2);
stem (a),xlabel('n'), ylabel('r(n)')
grid on
title('Sawtooth sequence');

7
OUTPUT:
Unit impulse function:
enter value of n for Impulse Function: 5

Unit Step function:

enter value of n for step function: 5

8
Unit Ramp function:

enter value of n for ramp function: 5

Exponential function:

enter value of n for exponential function: 5

9
Parabolic function:
enter value of n for parabola function: 5

Sinusoidal function:

10
Square function:
enter time of plot for square function: 20

Rectangle function:

11
Enter duty cycle for rectangle function: 75
Enter amplitude for rectangle function: 5

Triangular function:
Enter amplitude of the Triangular waveform: 5

Sawtooth function:

12
Enter amplitude of the sawtooth waveform: 5

RESULT:

13
LINEAR CONVOLUTION
AIM: To write MATLAB programs to find the linear convolution using MATLAB command, multiplication
method and matrix method.

SOFTWARE REQUIRED: MATLAB

THEORY:
The mathematical definition of convolution in discrete time domain is

y(n) = x(n) ∗ h(n) = ∑ 𝑥(𝑘)ℎ(𝑛 − 𝑘)


𝑘=−∞
Where x(n) is the input signal, h(n) is the impulse response, and y(n) is the output signal, ‘*’ denotes convolution.
If the input x(n) has N1 samples and h(n) has N2 samples then the output sequence y(n) will be a finite duration
sequence consisting of [N1+N2-1] samples. Here we multiply the terms of x(k) by the terms of the time shifted
h(n) and add them up.
In this equation x(k), h(n-k) and y(n) represents the input to the output from the system at the time ‘n’.
Here one of the inputs is shifted in time by a value every time it is multiplied with the other input signal. Linear
Convolution is quite often used as a method of implementing filters of various types.

ALORITHM:
1) Using MATLAB function
Step1: Enter the input sequence 𝑥(𝑛).
Step2: Enter the impulse response sequence ℎ(𝑛).
Step3: Perform the linear convolution between 𝑥(𝑘) and ℎ(𝑘) and obtain 𝑦(𝑛) using MATLAB
command CONV.
Step4: Plot x(𝑛), ℎ(𝑛) and 𝑦(𝑛).
2) Using Multiplication Method
Step1: Enter the input sequence 𝑥(𝑛).
Step2: Enter the impulse response sequence ℎ(𝑛).
Step3: Perform the zero padding such that the sequence length will be n+m-1, where n is length of
𝑥(𝑛), m is length of ℎ(𝑛).
Step4: Perform the linear convolution between 𝑥(𝑘) and ℎ(𝑘) using multiplication method and obtain
𝑦(𝑛).
Step5: Plot x(𝑛), ℎ(𝑛) and 𝑦(𝑛).
3) Using Matrix Method
Step1: Enter the input sequence 𝑥(𝑛).
Step2: Enter the impulse response sequence ℎ(𝑛).

14
Step3: Perform the zero padding such that the sequence length will be n+m-1, where n is length of
𝑥(𝑛), m is length of ℎ(𝑛).
Step4: Flip the one of the sequences. Let here ℎ(𝑛) sequence be flipped.
Step5: Obtain a square matrix by rotating the elements in the sequence and by row concatenating them.
Step6: Now transpose the matrix and multiply it with input sequence to obtain the convolved sequence
𝑦(𝑛).
Step7: Plot x(𝑛), ℎ(𝑛) and 𝑦(𝑛).

PROCEDURE:
1. Start the MATLAB program.
2. Create a NEW M-File from the file Menu.
3. An editor window open, type the program in it.
4. Save the file 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.

SOURCE CODE:
Using MATLAB Command:
%% Convolution using MATLAB Command (conv ())
clc
clear all
close all
%% Input sequence 1
x = input('enter values of 1st sequence: ');
n1 = length(x); % To get the length of the sequence
b = input('initial time of 1st signal: ');
nx= b:n1-1+b;
%% Input sequence 2
h = input('enter values of 2nd sequence: ');
n2 = length(h); % To get the length of the sequence
c = input('initial time of 2nd signal: ');
nh= c:n2-1+c;
%% Convolution
y= conv(x,h)
%% Plotting the sequences
% Sequence1
subplot(2,2,1);
stem(nx,x),xlabel('n'),ylabel('x(n)')
title('x(n)')
grid on
% Sequence2

15
subplot(2,2,2);
stem(nh,h),xlabel('n'),ylabel('h(n)')
title('h(n)')
grid on
% convolved signal
ny = b+c:1:n1+n2-2+b+c;
subplot(2,1,2);
stem(ny,y),xlabel('n'),ylabel('y(n)');
title('Convolution of sequences x(n) & h(n)');
grid on
Using Multiplication method:
%% Convolution using Multiplication method
clc
clear all
close all
%% Input sequence 1
x = input('enter values of 1st sequence: ');
n1 = length(x); % To get the length of the sequence
b = input('initial time of 1st signal: ');
nx= b:n1-1+b;
%% Input sequence 2
h = input('enter values of 2nd sequence: ');
n2 = length(h); % To get the length of the sequence
c = input('initial time of 2nd signal: ');
nh= c:n2-1+c;
%% Convolution
X = [x,zeros(1,n2)];
H = [h,zeros(1,n1)];
for i = 1:n1+n2-1
y(i) = 0;
for j = 1:n1
if (i-j+1>0)
y(i)=y(i)+X(j)*H(i-j+1);
end
end
end
%% Ploting the sequences
% Sequence1
subplot(2,2,1);
stem(nx,x),xlabel('n'),ylabel('x(n)')
title('x(n)')
grid on
% Sequence2
subplot(2,2,2);
stem(nh,h),xlabel('n'),ylabel('h(n)')

16
title('h(n)')
grid on
% convolved signal
ny = b+c:1:n1+n2-2+b+c;
subplot(2,1,2);
stem(ny,y),xlabel('n'),ylabel('y(n)');
title('Convolution of sequences x(n) & h(n)');
grid on
Using Matrix method:
%% Convolution using Matrix method
clc
clear all
close all
%% Input sequence 1
x = input('enter values of 1st sequence: ');
n1 = length(x);
b = input('initial time of 1st signal: ');
nx= b:n1-1+b;
%% Input sequence 2
h = input('enter values of 2nd sequence: ');
n2 = length(h);
c = input('initial time of 2nd signal: ');
nh= c:n2-1+c;
%% Convolution
X = [x,zeros(1,n2)];
H = [h,zeros(1,n1)];
H2 = fliplr(H);
for n3 = 1:n1+n2-1
H2 = [H2(end),H2(1:end-1)];
y(n3) = X * H2';
end
%% Ploting the sequences
% Sequence1
subplot(2,2,1);
stem(nx,x),xlabel('n'),ylabel('x(n)')
title('x(n)')
grid on
% Sequence2
subplot(2,2,2);
stem(nh,h),xlabel('n'),ylabel('h(n)')
title('h(n)')
grid on
% Convolved signal
ny = b+c:1:n1+n2-2+b+c;
subplot(2,1,2);

17
stem(ny,y),xlabel('n'),ylabel('y(n)');

title('Convolution of sequences x(n) & h(n)');


grid on

18
CALCULATIONS:

19
OUTPUT:
enter values of 1st sequence: [3,2,1,2]
initial time of 1st signal: 0
enter values of 2nd sequence: [1,2,1,2]
initial time of 2nd signal: -1

y=

3 8 8 12 9 4 4

RESULT:

20
CIRCULAR CONVOLUTION
AIM: To write MATLAB programs to find the circular convolution using MATLAB command, FFT method,
multiplication method and matrix method.

SOFTWARE REQUIRED: MATLAB

THEORY:
The convolution of two periodic sequences with period N is called circular convolution of two signals.
The circular convolution of x(n) and h(n) is denoted by
𝑁−1

𝑦(𝑛) = x(n)⨀h(n) = ∑ 𝑥(𝑘)ℎ((𝑛 − 𝑚))𝑁


𝑚=0
where ℎ((𝑛 − 𝑚))𝑁 is the reflected and circulated translated version of h(n).
Circular convolution in terms of DFT and IDFT is given as follows
𝑦(𝑛) = x(n) ⨀ h(n) = 𝐼𝐷𝐹𝑇 (𝐷𝐹𝑇(𝑥 (𝑛)). 𝐷𝐹𝑇(ℎ(𝑛)))
It can be performed only if both the sequences consist of equal number of samples. If the sequences are different
in length then convert the smaller size sequence to that of larger size by appending zeros.

ALORITHM:
1) Using MATLAB function
Step-1: Enter the input sequence 𝑥(𝑛).
Step-2: Enter the impulse response sequence ℎ(𝑛).
Step-3: Perform the circular convolution between 𝑥(𝑘) and ℎ(𝑘) and obtain 𝑦(𝑛) using MATLAB
command CCONV.
Step-4: Plot x(𝑛), ℎ(𝑛) and 𝑦(𝑛).
2) Using FFT and IDFT
Step-1: Enter the input sequence 𝑥(𝑛).
Step-2: Enter the impulse response sequence ℎ(𝑛).
Step-3: Perform FFT operation on 𝑥(𝑛) and ℎ(𝑛) to obtain 𝑋 (𝑘) and 𝐻 (𝑘) respectively using FFT
command.
Step-4: Multiply 𝑋(𝑘) and 𝐻 (𝑘) to obtain 𝑌(𝑘).
Step-5: Perform inverse FFT on 𝑌(𝑘) to obtain 𝑦(𝑛) using IFFT command
Step-6: Plot x(𝑛), ℎ(𝑛) and 𝑦(𝑛).
3) Using Multiplication \ Concentric circle Method
Step-1: Enter the input sequence 𝑥(𝑛).
Step-2: Enter the impulse response sequence ℎ(𝑛).

21
Step-3: Perform zero padding such that both the sequence length to ′𝑙′ where 𝑙 is maximum of length of
the sequences 𝑥(𝑛) and ℎ(𝑛).
Step-4: Perform the circular convolution between 𝑥(𝑛) and ℎ(𝑛) and obtain 𝑦(𝑛) using Multiplication \
Concentric circle method.
Step-5: Plot 𝑥(𝑛), ℎ(𝑛) and 𝑦(𝑛).
4) Using Matrix Method
Step-1: Enter the input sequence 𝑥(𝑛).
Step-2: Enter the impulse response sequence ℎ(𝑛).
Step-3: Perform zero padding such that both the sequence length to ′𝑙′ where 𝑙 is maximum of length of
the sequences 𝑥(𝑛) and ℎ(𝑛).
Step4: Flip the one of the sequences. Let here ℎ(𝑛) sequence be flipped.
Step5: Obtain a square matrix by rotating the elements in the sequence and by row concatenating them.
Step6: Now transpose the matrix and multiply it with input sequence to obtain the convolved sequence
𝑦(𝑛).
Step7: Plot x(𝑛), ℎ(𝑛) and 𝑦(𝑛).

PROCEDURE:
1. Start the MATLAB program.
2. Create a NEW M-File from the file Menu.
3. An editor window open, type the program in it.
4. Save the file 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.

SOURCE CODE:
Using MATLAB Command:
%% Circular Convolution using MATLAB Command (cconv ())
clc
clear all
close all
%% Input sequence 1
x = input('enter values of 1st sequence: ');
n1 = length(x);
nx= 0:n1-1;
%% Input sequence 2
h = input('enter values of 2nd sequence: ');
n2 = length(h);
nh= 0:n2-1;

22
%% Circular Convolution
l = max(n1,n2);
y = cconv(x,h,l)
%% Ploting the sequences
%Sequence1
subplot(2,2,1);
stem(nx,x),xlabel('n'),ylabel('x(n)')
title('x(n)')
grid on
%Sequence2
subplot(2,2,2);
stem(nh,h),xlabel('n'),ylabel('h(n)')
title('h(n)')
grid on
%convolved signal
ny = 0:l-1;
subplot(2,1,2);
stem(ny,y),xlabel('n'),ylabel('y(n)');
title('Circular convolution of sequences x(n) & h(n)');
grid on
Using FFT method:
%% Circular Convolution using FFT method
clc
clear all
close all
%% Input sequence 1
x = input('Enter values of 1st sequence: ');
n1 = length(x);
nx = 0:n1-1;
%% Input sequence 2
h = input('Enter values of 2nd sequence: ');
n2 = length(h);
nh = 0:n2-1;
%% Circular Convolution
y = ifft(fft(x) .* fft(h))
%% Ploting the sequences
% Sequence1
subplot(2,2,1);
stem(nx,x),xlabel('n'),ylabel('x(n)')
title('x(n)')
grid on
% Sequence2
subplot(2,2,2);
stem(nh,h),xlabel('n'),ylabel('h(n)')
title('h(n)')

23
grid on
% convolved signal
ny = 0:l-1;
subplot(2,1,2);
stem(ny,y),xlabel('n'),ylabel('y(n)');
title('Circular convolution of sequences x(n) & h(n)');
grid on
Using Concentric circle method:
%% Circular Convolution using Concentric Circle method
clc
clear all
close all
%% Input sequence 1
x = input('enter values of 1st sequence: ');
n1 = length(x);
nx= 0:n1-1;
%% Input sequence 2
h = input('enter values of 2nd sequence: ');
n2 = length(h);
nh= 0:n2-1;
%% Circular Convolution
l = max (n1,n2);
n3=n1-n2;
if(n3>0)
h=[h,zeros(1,N3)];
else
if (n3<0)
x=[x,zeros(1,-N3)];
else
end
end
for n=1:l;
y(n)=0;
for i=1:l;
j=n-i+1;
if(j<=0)
j=l+j;
end
y(n)=[y(n)+(x(i)*h(j))]
end
end
%% Ploting the sequences
% Sequence1
subplot(2,2,1);
stem(nx,x),xlabel('n'),ylabel('x(n)')

24
title('x(n)')
grid on
% Sequence2
subplot(2,2,2);
stem(nh,h),xlabel('n'),ylabel('h(n)')
title('h(n)')
grid on
% convolved signal
ny = 0:l-1;
subplot(2,1,2);
stem(ny,y),xlabel('n'),ylabel('y(n)');
title('Circular convolution of sequences x(n) & h(n)');
grid on
Using Matrix multiplication method:
%% Circular Convolution using Matrix multiplication
clc
clear all
close all
%% Input sequence 1
x = input('enter values of 1st sequence: ');
n1 = length(x);
nx= 0:n1-1;
%% Input sequence 2
h = input('enter values of 2nd sequence: ');
n2 = length(h);
nh= 0:n2-1;
%% Convolution
l = max (n1,n2);
n3=n1-n2;
if(n3>0)
h=[h,zeros(1,N3)];
else
if (n3<0)
x=[x,zeros(1,-N3)];
else
end
end
H = fliplr(h);
for n3 = 1:l
H = [H(end),H(1:end-1)];
y(n3) = x * H';
end
%% Ploting the sequences
% Sequence1
subplot(2,2,1);

25
stem(nx,x),xlabel('n'),ylabel('x(n)')
title('x(n)')
grid on
% Sequence2
subplot(2,2,2);
stem(nh,h),xlabel('n'),ylabel('h(n)')
title('h(n)')
grid on
% convolved signal
ny = 0:l-1;
subplot(2,1,2);
stem(ny,y),xlabel('n'),ylabel('y(n)');
title('Circular convolution of sequences x(n) & h(n)');
grid on

26
CALCULATIONS:

27
OUTPUT:
enter values of 1st sequence: [9,1,3,2]
enter values of 2nd sequence: [5,2,1,0,2]

y=

49 29 30 17 25

RESULT:

28
DISCRETE FOURIER TRANSFORM And INVERSE DISCRETE
FOURIER TRANSFORM
AIM: To write a MTLAB program for computation of DFT and IDFT using direct formula and FFT method.
SOFTWARE REQUIRED: MATLAB

THEORY:
DFT:
Discrete Fourier Transform (DFT) is used for performing frequency analysis of discrete time signal.
DFT gives a discrete frequency domain representation whereas the other transforms are continuous infrequency
domain. The N point DFT of discrete time signal x(n) is given by the equation
𝑁−1

𝑋 (𝑘 ) = ∑ 𝑥(𝑛)𝑒 −𝑗2𝜋𝑘𝑛/𝑁 ; 𝑘 = 0 𝑡𝑜 𝑁 − 1
𝑛=0
The inverse DFT allows us to recover the sequence x(n) from the frequency sample and it is given by
the equation
𝑁−1
1
𝑋(𝑛) = ∑ 𝑥(𝑘)𝑒 𝑗2𝜋𝑘𝑛/𝑁 ; 𝑛 = 0 𝑡𝑜 𝑁 − 1
𝑁
𝑘=0
FFT:
A Fast Fourier Transform (FFT) is an efficient algorithm to compute the DFT and IDFT. FFT are of
great importance to a wide variety of applications, from digital signal processing and solving partial differential
equations to algorithms for quick multiplication of large integers. Evaluating the sum of DFT directly wold take
O(N2) arithmetical operations. An FFT is an algorithm to compute the same result in only O(NlogN)
operations.

ALORITHM:
1) Using Direct Method
Step1: Get the input sequence.
Step2: Find the DFT of the input sequence using direct equation of DFT.
Step3: Find the IDFT of the sequence obtained in the above step using direct equation of IDFT.
Step4: Display the input sequence and output sequences obtained at each step using stem function.
2) Using FFT Method
Step1: Get the input sequence.
Step2: Find the DFT of the input sequence using FFT command.
Step3: Find the IDFT of the sequence using IFFT command.
Step4: Display the input sequence and output sequences obtained at each step using stem function.

PROCEDURE:
1. Start the MATLAB program.

29
2. Create a NEW M-File from the file Menu.
3. An editor window open, type the program in it.
4. Save the file 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.

SOURCE CODE:
Using direct formula:
%% Discrete Fourier Transform and Inverse Using Direct method
clc
clear all
close all
%% Sequence input
x = input('Enter the sequence : ')
%% DFT Operations
N = length(x);
for k = 1:N
X(k) = 0;
for n = 1:N
X(k) = X(k)+x(n)*exp(-1i*2*pi*(k-1)*(n-1)/N);
end
end
mag = abs(X);
pha = angle(X);
%% IDFT Operation
R = length(X);
for k = 1:N
x_idft(k) = 0;
for n = 1:N
x_idft(k) = x_idft(k)+X(n)*exp(i*2*pi*(k-1)*(n-1)/N);
end
end
x_idft = 1/R*(x_idft)
%% plotting the sequence
% Input signal
subplot(3,1,1);
stem(x),xlabel('n'),ylabel('x(n)');
title('Input signal');
grid on
% Magnitude Response
subplot(3,2,3);
stem(mag);xlabel('n'),ylabel('Magnitude');
title('Magnitude Response');
grid on

30
% Phase Response
subplot(3,2,4);
stem(pha);xlabel('n'),ylabel('phase');
title('phase Response');
grid on
% Inverse DFT
subplot(3,1,3);
stem(x_idft),xlabel('n'),ylabel('x(n)');
title('Inverse DFT');
grid on
Using FFT method:
%% Discrete Fourier Transform and Inverse Using Direct method
clc
clear all
close all
%% Sequence input
x = input('Enter the sequence : ')
%% DFT Operations
X = fft(x)
mag = abs(X);
pha = angle(X);
%% IDFT Operation
x_idft = ifft(X)
%% plotting the sequence
% Input signal
subplot(3,1,1);
stem(x),xlabel('n'),ylabel('x(n)');
title('Input signal');
grid on
% Magnitude Response
subplot(3,2,3);
stem(mag);xlabel('n'),ylabel('Magnitude');
title('Magnitude Response');
grid on
% Phase Response
subplot(3,2,4);
stem(pha);xlabel('n'),ylabel('phase');
title('phase Response');
grid on
% Inverse DFT
subplot(3,1,3);
stem(x_idft),xlabel('n'),ylabel('x(n)');
title('Inverse DFT');
grid on

31
32
CALCULATIONS:

33
OUTPUT:
Enter the sequence : [1 2 3 4 4 3 2 1]
x=

1 2 3 4 4 3 2 1

X=

20.0000 + 0.0000i -5.8284 - 2.4142i 0.0000 + 0.0000i -0.1716 - 0.4142i 0.0000 +


0.0000i -0.1716 + 0.4142i 0.0000 + 0.0000i -5.8284 + 2.4142i

x_idft =
1 2 3 4 4 3 2 1

34
RESULT:

35
ADDITION OF SINUSOIDAL SIGNALS
AIM: To write MATLAB programs to perform addition of given sinusoidal signals and to perform one of the
applications i.e. Amplitude Modulation.

SOFTWARE REQUIRED: MATLAB

THEORY:
Addition of signals is nothing but addition of amplitude of two or more signals at each instance of time or any
other independent variable which are common between the signals.
Amplitude modulation (AM) is one of the conventional modulation techniques to transmit signals using
a carrier wave. The amplitude or the strength of a high frequency carrier wave is changed in accordance with the
amplitude of message signal. Let’s consider carrier wave as Ac*sin(2пfmt) and message signals Am*sin(2пfmt)
where Ac and Am are amplitudes of carrier and message signals respectively. Then the modulated signal will be
(Ac+Am*sin(2пfmt))*sin(2п fct).
Modulation Index or the Modulation Depth is the one of the most common term that used along with
modulation techniques. Here in A.M., it is the measure of amplitude variation surrounding an unmodulated
carrier. It is the ratio of the amplitude of message signal to the amplitude of the carrier signal. In terms of
modulation index (m = Am/Ac) the equation of the modulated signal becomes (1+m* sin(2пfmt)) * Acsin(2п fct).

ALORITHM:
1) Sum of Sinusoidal
Step1: Enter the time period, frequency to generate the sinusoidal signals.
Step2: Perform addition operation on the given signals.
Step3: Plot the graphs.
2) Application (Amplitude Modulation)
Step1: Enter the time period.
Step2: Generate the message signal by taking the frequency as input.
Step3: Generate the carrier signal by taking the frequency as input.
Step4: Perform amplitude modulation by using 𝑆 = (1 + 𝑢 ∗ 𝑚) .∗ 𝑐
where, s − modulated signal
u − modulation index
m − message signal
c – carrier signal

Step5: Plot the graphs.

PROCEDURE:
1. Start the MATLAB program.
2. Create a NEW M-File from the file Menu.

36
3. An editor window open, type the program in it.
4. Save the file 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.

SOURCE CODE:
Addition of signals:
%% Addition of Sinusoidal signals
clc
clear all
close all
t=input('Enter time period for signals: ');
t1=linspace(0,t,1000);
%% Signal 1
fm=input('Enter Signal-1 frequency: ');
Am=input('Enter signal-1 amplitude: ');
m = Am.*sin(2*pi*fm*t1);
%% Signal 2
fc=input('Enter Signal-2 frequency: ');
Ac=input('Enter signal-2 amplitude: ');
c = Ac.*sin(2*pi*fc*t1);
%% Addition of signals
s = m+c;
%% Plotting the signals
% Signal 1
figure;
subplot(3,1,1);
plot(t1,m),xlabel('time'),ylabel('amplitude');
title('Signal 1');
% Signal 2
subplot(3,1,2);
plot(t1,c),xlabel('time'),ylabel('amplitude');
title('Signal 2');
% Added signal
subplot(3,1,3)
plot(s),xlabel('time'),ylabel('amplitude');
title('Summed signal');
Application of addition of signals:
%% Application of Addition of Sinusoidal signals
clc
clear all
close all
t=input('Enter time period for signals: ');

37
t1=linspace(0,t,1000);
%% Message Signal
fm=input('Enter message signal frequency: ');
m = sin(2*pi*fm*t1);
%% Carrier Signal
fc=input('Enter carrier signal frequency: ');
c = sin(2*pi*fc*t1);
%% Modulation
u = input('Enter modulation index: ');
s = (1+(u.*m)).*c;
%% Plotting the signals
% message signal
subplot(2,1,1);
plot(t1,m),xlabel('t'),ylabel('amplitude');
title('Message signal');
% Carrier Signal
subplot(2,1,2);
plot(t1,c),xlabel('t'),ylabel('amplitude');
title('Carrier Signal');
% modulated signal
figure;
plot(t1,s),xlabel('t'),ylabel('amplitude');
title('Modulated signal');

38
OUTPUT:
Output for sum of sinusoidal:
enter time period for signals: 10
enter Signal-1 frequency: 500
enter signal-1 amplitude: 10
enter Signal-2 frequency: 5000
enter signal-2 amplitude: 5

39
Output for Application of sum of sinusoidal:
enter time period for signals: 10
enter message signal frequency: 500
enter carrier signal frequency: 5000

enter modulation index: 1

40
enter modulation index: 0.7

enter modulation index: 1.7

41
RESULT:

42
CYCLE - 2

43
DESIGN OF ANALOG BUTTERWORTH FILTER
AIM: To design and implement the Analog Butterworth filters i.e. Low pass and High pass filter for the given
specifications using MATLAB.

SOFTWARE REQUIRED: MATLAB

THEORY:
Butterworth filters are causal in nature and of various orders, the lowest order being the best (shortest) in
the time domain, and the higher orders being better in the frequency domain. Butterworth or maximally flat filters
have a monotonic amplitude frequency response which is maximally flat at zero frequency response and the
amplitude frequency response decreases logarithmically with increasing frequency. A Butterworth filter is
characterized by its magnitude frequency response,
1
|𝐻 (𝑗Ω)| = 1/2
Ω 2𝑁
[1 + (Ω ) ]
𝑐

Where N is the order of the filter and Ωc is defined as the cut-off frequency where the filter magnitude is 1/√2
times the dc gain (Ω=0).
As N ∞, the low pass filter is said to be normalized. All types of filters namely – Low pass, High
pass, Band pass and Band elimination filters can be derived from the Normalized Low pass filter.
Steps to design an analog Butterworth low pass filter:
1. From the given specification find the order of the filter N.
2. Round off it to the next highest integer.
3. Find the transfer function H(s) for Ωc = 1 rad/sec for the value of N.
4. Calculate the value of cut-off frequency Ωc.
𝑆
5. Find the transfer function Ha(s) for the above value of Ωc by substituting 𝑆 → in H(s).
Ωc

ALORITHM:
Step-1: Enter the passband ripple (𝑟𝑝 ) and stopband ripple (𝑟𝑠 ).

Step-2: Enter the passband frequency (𝑓𝑝 ) and stopband frequency (𝑓𝑠 ).

For LPF −> 𝑓𝑝 < 𝑓𝑠 and for 𝐻𝑃𝐹 −> 𝑓𝑝 > 𝑓𝑠 .

Step-3: Get the sampling frequency (𝑓).


Step-4: Calculate the normalized passband, stopband frequencies i.e. 𝑤𝑝 , 𝑤𝑠 respectively.

𝑤𝑝 = 2𝑓𝑝 ⁄𝑓 ; 𝑤𝑠 = 2𝑓𝑠 ⁄𝑓.

Step-5: Calculate the order and 3dB cut-off frequency using


[𝑛, 𝑤𝑛 ] = 𝑏𝑢𝑡𝑡𝑜𝑟𝑑(𝑤𝑝 , 𝑤𝑠 , 𝑟𝑝 , 𝑟𝑠, ′s′ )

Step-6: Find the window coefficients for low pass filter using

44
[𝑏, 𝑎] = butter (n, 𝑤𝑛 , ‘low’, ‘s’)

(𝑜𝑟)
Find the window coefficients for high pass filter using
[𝑏, 𝑎] = butter (n, 𝑤𝑛 , ‘high’, ‘s’)

Step-7: Find the frequency response using ‘freqs( )’ function.


Step-8: Calculate and plot the magnitude response and phase response.

PROCEDURE:
1. Start the MATLAB program.
2. Create a NEW M-File from the file Menu.
3. An editor window open, type the program in it.
4. Save the file 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.

SOURCE CODE:
Low Pass Filter:
%% BUTTERWORT LOW PASS FILTER DESIGN
clc;
clear all;
close all;
%% Filter Specifications
disp('Enter Filter Specifications');
f = input('Enter Sampling frequency in Hz : ');
% inputs for stop band
fs = input('Stopband frequency in Hz : ');
ws = 2*fs/f;
rs = input('Stopband attenuation in dB : ');
% inputs for pass band
fp = input('Passband frequency in Hz : ');
wp = 2*fp/f;
rp = input('Passband attenuation in dB : ');
%% Order and cutoff frequency of filter
[n,wn] = buttord(wp,ws,rp,rs,'s');
%% Response of filter
[b,a] = butter(n,wn,'s')
%% Magnitude and Phase response
w = 0:0.01:pi;
[h1,om] = freqs(b,a,w);

45
m = 20*log10(abs(h1));
an = angle(h1);
%% Plotting
% MAGNITUDE PLOT
subplot(2,1,1);
plot(om/pi,m),xlabel('NORMALISED FREQUENCY'),ylabel('GAIN IN
DB');
title('MAGNITUDE RESPONSE OF BUTTERWORTH LOW PASS FILTER');
grid on;
% PHASE PLOT
subplot(2,1,2);
plot(om/pi,an),xlabel('NORMALISED FREQUENCY'),ylabel('PHASE');
title('PHASE RESPONSE OF BUTTERWORTH LOW PASS FILTER');
grid on;

High Pass Filter:


%% BUTTERWORT HIGH PASS FILTER DESIGN
clc;
clear all;
close all;
%% Filter Specifications
disp('Enter Filter Specifications');
f = input('Enter Sampling frequency in Hz : ');
% inputs for stop band
fs = input('Stopband frequency in Hz : ');
ws = 2*fs/f;
rs = input('Stopband attenuation in dB : ');
% inputs for pass band
fp = input('Passband frequency in Hz : ');
wp = 2*fp/f;
rp = input('Passband attenuation in dB : ');
%% Order and cutoff frequency of filter
[n,wn] = buttord(wp,ws,rp,rs,'s');
%% Response of filter
[b,a] = butter(n,wn,'high','s')
%% Magnitude and Phase response
w = 0:0.01:pi;
[h1,om] = freqs(b,a,w);
m = 20*log10(abs(h1));
an = angle(h1);
%% Plotting
% MAGNITUDE PLOT
subplot(2,1,1);
plot(om/pi,m),xlabel('NORMALISED FREQUENCY'),ylabel('GAIN IN
DB');

46
title('MAGNITUDE RESPONSE OF BUTTERWORTH HIGH PASS FILTER');
grid on;
% PHASE PLOT
subplot(2,1,2);
plot(om/pi,an),xlabel('NORMALISED FREQUENCY'),ylabel('PHASE');
title('PHASE RESPONSE OF BUTTERWORTH HIGH PASS FILTER');
grid on;

47
CALCULATIONS:

48
OUTPUT:
Output for LPF:
Enter Filter Specifications
Enter Sampling frequency in Hz : 200
Stopband frequency in Hz : 3000
Stopband attenuation in dB : 10
Passband frequency in Hz : 2000
Passband attenuation in dB : 2

b=
1.0e+05 *
0 0 0 0 2.7000
a=
1.0e+05 *

49
0.0000 0.0006 0.0177 0.3095 2.7000

50
Output for HPF:
Enter Filter Specifications
Enter Sampling frequency in Hz : 200
Stopband frequency in Hz : 2000
Stopband attenuation in dB : 10
Passband frequency in Hz : 3000
Passband attenuation in dB : 2

b=
1 0 0 0 0
a=
1.0e+05 *
0.0000 0.0007 0.0237 0.4765 4.8000

RESULT:

51
DESIGN OF ANALOG CHEBYSHEV FILTER
AIM: To design and implement the Analog Chebyshev filters i.e. type-1 Low pass\High pass filter and type-2
Low pass\High pass filter for the given specifications using MATLAB.

SOFTWARE REQUIRED: MATLAB

THEORY:
Chebyshev filters are equi-ripple in either the passband or stopband. Hence the magnitude response
oscillates between the permitted minimum and maximum values in the band a number of times depending upon
the order of filters. There are two types of Chebyshev filters. The Chebyshev I filter is equi-ripple in passband
and monotonic in the stopband, whereas Chebyshev II is just the opposite. The Chebyshev low-pass filter has a
magnitude response given by
−1
2 2 2 Ω
( )
|𝐻 𝑗Ω | = [1 + 𝜀 𝐶𝑁 ( )]
Ω𝑐
where 𝝴 is a parameter related to the ripple present in the passband
CN(x) is given by
cos(𝑁 cos −1 𝑥 ) , |𝑥| ≤ 1, 𝑃𝑎𝑠𝑠𝑏𝑎𝑛𝑑
𝐶𝑁 (𝑥 ) = {
cos(𝑁 cosh−1 𝑥 ) , |𝑥| > 1, 𝑆𝑡𝑜𝑝𝑏𝑎𝑛𝑑
The magnitude response has equi-ripple pass band and maximally flat stop band. By increasing the filter
order N, the Chebyshev response approximates the ideal response. The phase response of the Chebyshev filter is
more non-linear than the Butter worth filter for a given filter length N.
Steps to design an analog Chebyshev type-1 low pass filter:
1. From the given specifications find the order of filer N.
2. Round-off it to the next higher integer.
3. Using the following formula find the values of ‘a’ and ‘b’, which are minor and major axis of ellipse
respectively.
[𝜇 1⁄𝑁 −𝜇 −1⁄𝑁 ] [𝜇 1⁄𝑁 +𝜇 −1⁄𝑁 ]
𝑎 = Ω𝑝 ; 𝑏 = Ω𝑝 .
2 2
where 𝜇 = 𝜖 −1 + √𝜖 −2 + 1 ;
𝜖 = √100.1𝛼𝑝 − 1 ;
Ωp is passband frequency;
αp is maximum allowable attenuation in the passband.
[For normalized Chebyshev filter Ωp = 1rad/sec]
4. Calculate the poles of Chebyshev filter which lies on an ellipse by using formula
𝑆ₖ = 𝑎 cos ∅ₖ + 𝑏 sin ∅ₖ , 𝑘 = 1 𝑡𝑜 𝑁
𝜋 2𝑘−1
where ∅ₖ = 2 + 2𝑁 𝜋 , 𝑘 = 1 𝑡𝑜 𝑁
5. Find the denominator polynomial of the transfer function using the above poles.
6. The numerator of the transfer function depends on the value of N.
a. For ‘N’ odd substitute s = 0 in the denominator polynomial and find the value. This is equal to
the numerator of the transfer function.
b. For ‘N’ even substitute s = 0 in the denominator polynomial and divide the result by √1 + 𝜖 2 .
This is equal to the numerator.

52
ALORITHM:
Step-1: Enter the passband ripple (𝑟𝑝 ) and stopband ripple (𝑟𝑠 ).

Step-2: Enter the passband frequency (𝑓𝑝 ) and stopband frequency (𝑓𝑠 ).

For LPF −> 𝑓𝑝 < 𝑓𝑠 and for 𝐻𝑃𝐹 −> 𝑓𝑝 > 𝑓𝑠 .

Step-3: Get the sampling frequency (𝑓).


Step-4: Calculate the normalized passband, stopband frequencies i.e. 𝑤𝑝 , 𝑤𝑠 respectively.

𝑤𝑝 = 2𝑓𝑝 ⁄𝑓 ; 𝑤𝑠 = 2𝑓𝑠 ⁄𝑓.

Step-5: Calculate the order and 3dB cut-off frequency Chebyshev-I LPF and HPF using [𝑛, 𝑤𝑛 ] =
𝑏𝑢𝑡𝑡𝑜𝑟𝑑(𝑤𝑝 , 𝑤𝑠 , 𝑟𝑝 , 𝑟𝑠, ′s′ )

(or)
Calculate the order and 3dB cut-off frequency for Chebyshev-II LPF and HPF using
[𝑛, 𝑤𝑛 ] = 𝑐ℎ𝑒𝑏2𝑜𝑟𝑑(𝑤𝑝 , 𝑤𝑠 , 𝑟𝑝 , 𝑟𝑠, ′s′ )

Step-6: Find the window co-efficient for Chebyshev-I of low pass and high pass using
[𝑏, 𝑎] = cheby1 (n, 𝑟𝑝 , 𝑤𝑛 , ‘low’, ‘s’) and [𝑏, 𝑎] = cheby1 (n, 𝑟𝑝 , 𝑤𝑛 , ‘high’, ‘s’) respectively

(or)
Find the window co-efficient for Chebyshev-II of low pass and high pass using
[𝑏, 𝑎] = cheby2 (n, 𝑟𝑠 , 𝑤𝑛 , ‘low’, ‘s’) and [𝑏, 𝑎] = cheby2 (n, 𝑟𝑠 , 𝑤𝑛 , ‘high’, ‘s’) respectively.
Step-7: Find the frequency response using ‘freqs( )’ function.
Step-8: Calculate and plot the magnitude and phase response.

PROCEDURE:
1. Start the MATLAB program.
2. Create a NEW M-File from the file Menu.
3. An editor window open, type the program in it.
4. Save the file 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.

SOURCE CODE:
Type-1 Low Pass Filter:
%% Chebyshev type-1 lowpass filter
clc;
clear all;
close all;

53
%% Filter Specifications
disp('Enter Filter Specifications');
% Sampling Frequency
f = input('Enter Sampling frequency in Hz : ');
% Inputs for pass band
fp= input('Passband frequency in Hz : ');
wp = 2*fp/f;
rp= input('Passband attenuation in dB : ');
% Inputs for stop band
fs= input('Stopband frequency in Hz : ');
ws = 2*fs/f;
rs= input('Stopband attenuation in dB : ');
%% Order and cutoff frequency of filter
[N,WN]=cheb1ord(wp,ws,rp,rs,'s');
%% Response of filter
[B,A]=cheby1(N,rp,WN,'low','s')
%% Magnitude and Phase response
w=0:0.01:pi;
[h,om]=freqs(B,A,w);
m=20*log10(abs(h));
an=angle(h);
%% Plotting
% MAGNITUDE PLOT
subplot(2,1,1);
plot(om/pi,m),xlabel('NORMALISED FREQUENCY'),ylabel('GAIN IN
DB');
title('MAGNITUDE RESPONSE OF CHEBYSHEV TYPE-1 LOW PASS FILTER');
grid on;
% PHASE PLOT
subplot(2,1,2);
plot(om/pi,an),xlabel('NORMALISED FREQUENCY'),ylabel('PHASE');
title('PHASE RESPONSE OF CHEBYSHEV TYPE-1 LOW PASS FILTER');
grid on;

Type-1 High Pass Filter:


%% Chebyshev type-1 highpass filter
clc;
clear all;
close all;
%% Filter Specifications
disp('Enter Filter Specifications');
% Sampling Frequency
f = input('Enter Sampling frequency in Hz : ');
% Inputs for pass band
fp= input('Passband frequency in Hz : ');

54
wp = 2*fp/f;
rp= input('Passband attenuation in dB : ');
% Inputs for stop band
fs= input('Stopband frequency in Hz : ');
ws = 2*fs/f;
rs= input('Stopband attenuation in dB : ');
%% Order and cutoff frequency of filter
[N,WN]=cheb1ord(wp,ws,rp,rs,'s');
%% Response of filter
[B,A]=cheby1(N,rp,WN,'high','s')
%% Magnitude and Phase response
w=0:0.01:pi;
[h,om]=freqs(B,A,w);
m=20*log10(abs(h));
an=angle(h);
%% Plotting
% MAGNITUDE PLOT
subplot(2,1,1);
plot(om/pi,m),xlabel('NORMALISED FREQUENCY'),ylabel('GAIN IN
DB');
title('MAGNITUDE RESPONSE OF CHEBYSHEV TYPE-1 HIGH PASS
FILTER');
grid on;
% PHASE PLOT
subplot(2,1,2);
plot(om/pi,an),xlabel('NORMALISED FREQUENCY'),ylabel('PHASE');
title('PHASE RESPONSE OF CHEBYSHEV TYPE-1 HIGH PASS FILTER');
grid on;

Type-2 Low Pass Filter:


%% Chebyshev type-2 lowpass filter
clc;
clear all;
close all;
%% Filter Specifications
disp('Enter Filter Specifications');
% Sampling Frequency
f = input('Enter Sampling frequency in Hz : ');
% Inputs for pass band
fp= input('Passband frequency in Hz : ');
wp = 2*fp/f;
rp= input('Passband attenuation in dB : ');
% Inputs for stop band
fs= input('Stopband frequency in Hz : ');
ws = 2*fs/f;

55
rs= input('Stopband attenuation in dB : ');
%% Order and cutoff frequency of filter
[N,WN]=cheb2ord(wp,ws,rp,rs,'s');
%% Response of filter
[B,A]=cheby2(N,rs,WN,'low','s')
%% Magnitude and Phase response
w=0:0.01:pi;
[h,om]=freqs(B,A,w);
m=20*log10(abs(h));
an=angle(h);
%% Plotting
% MAGNITUDE PLOT
subplot(2,1,1);
plot(om/pi,m),xlabel('NORMALISED FREQUENCY'),ylabel('GAIN IN
DB');
title('MAGNITUDE RESPONSE OF CHEBYSHEV TYPE-2 LOW PASS FILTER');
grid on;
% PHASE PLOT
subplot(2,1,2);
plot(om/pi,an),xlabel('NORMALISED FREQUENCY'),ylabel('PHASE');
title('PHASE RESPONSE OF CHEBYSHEV TYPE-2 LOW PASS FILTER');
grid on;

Type-2 High Pass Filter:


%% Chebyshev type-2 highpass filter
clc;
clear all;
close all;
%% Filter Specifications
disp('Enter Filter Specifications');
% Sampling Frequency
f = input('Enter Sampling frequency in Hz : ');
% Inputs for pass band
fp= input('Passband frequency in Hz : ');
wp = 2*fp/f;
rp= input('Passband attenuation in dB : ');
% Inputs for stop band
fs= input('Stopband frequency in Hz : ');
ws = 2*fs/f;
rs= input('Stopband attenuation in dB : ');
%% Order and cutoff frequency of filter
[N,WN]=cheb2ord(wp,ws,rp,rs,'s');
%% Response of filter
[B,A]=cheby2(N,rs,WN,'high','s')
%% Magnitude and Phase response

56
w=0:0.01:pi;
[h,om]=freqs(B,A,w);
m=20*log10(abs(h));
an=angle(h);
%% Plotting
% MAGNITUDE PLOT
subplot(2,1,1);
plot(om/pi,m),xlabel('NORMALISED FREQUENCY'),ylabel('GAIN IN
DB');
title('MAGNITUDE RESPONSE OF CHEBYSHEV TYPE-2 HIGH PASS
FILTER');
grid on;
% PHASE PLOT
subplot(2,1,2);
plot(om/pi,an),xlabel('NORMALISED FREQUENCY'),ylabel('PHASE');
title('PHASE RESPONSE OF CHEBYSHEV TYPE-2 HIGH PASS FILTER');
grid on;

57
CALCULATIONS:

58
OUTPUT:
Output of Type-1 LPF:
Enter Filter Specifications
Enter Sampling frequency in Hz : 7000
Passband frequency in Hz : 1500
Passband attenuation in dB : 0.15
Stopband frequency in Hz : 3000
Stopband attenuation in dB : 60
B=
1.0e-04 *
0 0 0 0 0 0 0 0 0.4743
A=
1.0000 0.6621 0.5865 0.2492 0.1026 0.0268 0.0056 0.0007 0.0000

59
Output of Type-1 HPF:
Enter Filter Specifications
Enter Sampling frequency in Hz : 7000
Passband frequency in Hz : 1500
Passband attenuation in dB : 0.15
Stopband frequency in Hz : 3000
Stopband attenuation in dB : 60
B=

0.9829 0 0 0 0 0 0 0 0
A=

1.0000 2.7347 3.9389 3.4404 2.4206 1.0795 0.4667 0.0968 0.0268

60
Output of Type-2 LPF:
Enter Filter Specifications
Enter Sampling frequency in Hz : 7000
Passband frequency in Hz : 1500
Passband attenuation in dB : 0.15
Stopband frequency in Hz : 3000
Stopband attenuation in dB : 60
B=

0.0010 -0.0000 0.0180 -0.0000 0.0507 -0.0000 0.0457 -0.0000 0.0129


A=

1.0000 2.8323 4.0110 3.6975 2.4286 1.1697 0.4086 0.0966 0.0129

Output of Type-2 HPF:

61
Enter Filter Specifications
Enter Sampling frequency in Hz : 7000
Passband frequency in Hz : 1500
Passband attenuation in dB : 0.15
Stopband frequency in Hz : 3000
Stopband attenuation in dB : 60
B=

1.0000 0.0000 1.1261 0.0000 0.3963 0.0000 0.0446 0.0000 0.0008


A=

1.0000 4.2295 10.0703 16.2309 18.9734 16.2640 9.9336 3.9494 0.7851

62
RESULT:

63
DESIGN OF IIR DIGITAL FILTERS
AIM: To design and implementation of IIR filter (Low pass\High pass) to meet the given specifications using
Bilinear transformation and Impulse Invariance technique.

SOFTWARE REQUIRED: MATLAB

THEORY:
The design of IIR filters is closely related to the design of analogue filters, which is a widely studied topic. An
analogue filter is usually designed and a transformation is carried out into the digital domain. IIR filters are
designed essentially by the Impulse Invariance or the Bilinear Transformation method.
Bilinear Transformation:
The Bilinear Transformation method overcomes the effect of aliasing that is caused due to the analogue
frequency response containing components at or beyond the Nyquist Frequency. The bilinear transform is a
method of compressing the infinite, straight analogue frequency axis to a finite one long enough to wrap around
the unit circle once only. This is also sometimes called frequency warping. This introduces a distortion in the
frequency. This is undone by pre-warping the critical frequencies of the analogue filter (cut-off frequency, centre
frequency) such that when the analogue filter is transformed into the digital filter, the designed digital filter will
meet the desired specifications.
Steps to design filter using bilinear transformation technique:
2 𝜔
1. From the given specification, find prewarping analogue frequencies using the formula Ω = 𝑇 tan 2 .
2. Using the analogue frequencies find H(s) of the analogue filter.
3. Select the sampling rate of the digital filter, call it ‘T’ seconds per sample.
2 1−𝑍 −1
4. Substitute 𝑠 = 𝑇 1+𝑍 −1 into the transfer function found in step-2
Impulse Invariance
In impulse invariance method the IIR filter is designed such hat unit impulse response h(n) of digital filter
is the sampled version of the impulse response of the analogue filter.
Steps to design filter using Impulse Invariance technique:
1. For the given specifications, find Ha(s), the transfer function of an analogue filter.
2. Select the sampling rate of the digital filter, T second per sample.
3. Express the analogue filter transfer function as the sum of single-pole filter.
𝑁
𝑐ₖ
𝐻ₐ(𝑠) = ∑
𝑠 − 𝑝ₖ
𝑘=1
4. Compute the Z-transform of the digital filter by using the formula
𝑁
𝑐ₖ
𝐻 (𝑧 ) = ∑ −1
1 − 𝑒 𝑝ₖ𝑇𝑧
𝑘=1

ALORITHM:
Step-1: Enter the passband ripple (𝑟𝑝 ) and stopband ripple (𝑟𝑠 ).

Step-2: Enter the passband frequency (𝑓𝑝 ) and stopband frequency (𝑓𝑠 ).

For LPF −> 𝑓𝑝 < 𝑓𝑠 and for 𝐻𝑃𝐹 −> 𝑓𝑝 > 𝑓𝑠 .

64
Step-3: Get the sampling frequency (𝑓).
Step-4: Calculate the normalized passband, stopband frequencies i.e. 𝑤𝑝 , 𝑤𝑠 respectively.

𝑤𝑝 = 2𝑓𝑝 ⁄𝑓;𝑤𝑠 = 2𝑓𝑠 ⁄𝑓

Step-5: Perform pre-warping operation.


Step-6: Calculate the order and 3dB cut-off frequency using appropriate formula for analog filter.
Step-7: Find the analog filter window coefficients using appropriate formulas (or) function.
Step-8: Now convert the analog filter into digital filters by using ‘bilinear ( )’ function for bilinear
transformation and ‘impinvar ( )’ function for impulse invariant method.
Step-9: Find the frequency response using ‘freqz( )’ function.
Step-10: Calculate and plot the magnitude and phase response.

PROCEDURE:
1. Start the MATLAB program.
2. Create a NEW M-File from the file Menu.
3. An editor window open, type the program in it.
4. Save the file 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.

SOURCE CODE:
Bilinear transformation Low pass filter:
%% Bilinear transformation (LPF)
clc;
clear all;
close all;
%% Filter Specifications
disp('Enter Filter Specifications');
f = input('Enter Sampling frequency in Hz : ');
% inputs for stop band
fs = input('Stopband frequency in Hz : ');
ws = 2*fs/f;
rs = input('Stopband attenuation in dB : ');
% inputs for pass band
fp = input('Passband frequency in Hz : ');
wp = 2*fp/f;
rp = input('Passband attenuation in dB : ');
%% Prewarping
t = 1/f;
omegap = (2/t)*tan(wp/2);

65
omegas = (2/t)*tan(ws/2);
%% Cutoff frequency and order of filter
[n,wn] = buttord(omegap, omegas,rp,rs,'s');
%% Response of filter
[c,d] = butter(n,wn,'low','s');
%% Bilinear Transformation
[b,a]= bilinear(c,d,f)
%% Magnitude and Phase response
w = 0:0.01:pi;
[h1,om] = freqz(b,a,w);
m = 20*log10(abs(h1));
an = angle(h1);
%% Plotting
% MAGNITUDE PLOT
subplot(2,1,1);
plot(om/pi,m),xlabel('NORMALISED FREQUENCY'),ylabel('GAIN IN
DB');
title('MAGNITUDE RESPONSE');
grid on;
% PHASE PLOT
subplot(2,1,2);
plot(om/pi,an),xlabel('NORMALISED FREQUENCY'),ylabel('PHASE');
title('PHASE RESPONSE');
grid on;

Bilinear transformation High pass filter:


%% Bilinear transformation (HPF)
clc;
clear all;
close all;
%% Filter Specifications
disp('Enter Filter Specifications');
f = input('Enter Sampling frequency in Hz : ');
% inputs for stop band
fs = input('Stopband frequency in Hz : ');
ws = 2*fs/f;
rs = input('Stopband attenuation in dB : ');
% inputs for pass band
fp = input('Passband frequency in Hz : ');
wp = 2*fp/f;
rp = input('Passband attenuation in dB : ');
%% Prewarping
t = 1/f;
omegap = (2/t)*tan(wp/2);
omegas = (2/t)*tan(ws/2);

66
%% Cutoff frequency and order of filter
[n,wn] = buttord(omegap, omegas,rp,rs,'s');
%% Response of filter
[c,d] = butter(n,wn,'high','s');
%% Bilinear Transformation
[b,a]= bilinear(c,d,f)
%% Magnitude and Phase response
w = 0:0.01:pi;
[h1,om] = freqz(b,a,w);
m = 20*log10(abs(h1));
an = angle(h1);
%% Plotting
% MAGNITUDE PLOT
subplot(2,1,1);
plot(om/pi,m),xlabel('NORMALISED FREQUENCY'),ylabel('GAIN IN
DB');
title('MAGNITUDE RESPONSE');
grid on;
% PHASE PLOT
subplot(2,1,2);
plot(om/pi,an),xlabel('NORMALISED FREQUENCY'),ylabel('PHASE');
title('PHASE RESPONSE');
grid on;

Impulse Invariance technique Low pass filter:


%% Impulse invariance method (LPF)
clc;
clear all;
close all;
%% Filter Specifications
disp('Enter Filter Specifications');
f = input('Enter Sampling frequency in Hz : ');
% inputs for stop band
fs = input('Stopband frequency in Hz : ');
ws = 2*fs/f;
rs = input('Stopband attenuation in dB : ');
% inputs for pass band
fp = input('Passband frequency in Hz : ');
wp = 2*fp/f;
rp = input('Passband attenuation in dB : ');
%% Prewarping
t = 1/f;
omegap = wp/t;
omegas = ws/t;
%% Cutoff frequency and order of filter

67
[n,wn] = buttord(omegap, omegas,rp,rs,'s');
%% Response of filter
[c,d] = butter(n,wn,'low','s');
%% Bilinear Transformation
[b,a]= impinvar(c,d,f)
%% Magnitude and Phase response
w = 0:0.01:pi;
[h1,om] = freqz(b,a,w);
m = 20*log10(abs(h1));
an = angle(h1);
%% Plotting
% MAGNITUDE PLOT
subplot(2,1,1);
plot(om/pi,m),xlabel('NORMALISED FREQUENCY'),ylabel('GAIN IN
DB');
title('MAGNITUDE RESPONSE');
grid on;
% PHASE PLOT
subplot(2,1,2);
plot(om/pi,an),xlabel('NORMALISED FREQUENCY'),ylabel('PHASE');
title('PHASE RESPONSE');
grid on;

Impulse Invariance technique High pass filter:


%% Impulse invariance method (HPF)
clc;
clear all;
close all;
%% Filter Specifications
disp('Enter Filter Specifications');
f = input('Enter Sampling frequency in Hz : ');
% inputs for stop band
fs = input('Stopband frequency in Hz : ');
ws = 2*fs/f;
rs = input('Stopband attenuation in dB : ');
% inputs for pass band
fp = input('Passband frequency in Hz : ');
wp = 2*fp/f;
rp = input('Passband attenuation in dB : ');
%% Prewarping
t = 1/f;
omegap = wp/t;
omegas = ws/t;
%% Cutoff frequency and order of filter
[n,wn] = buttord(omegap, omegas,rp,rs,'s');

68
%% Response of filter
[c,d] = butter(n,wn,'high','s');
%% Bilinear Transformation
[b,a]= impinvar(c,d,f)
%% Magnitude and Phase response
w = 0:0.01:pi;
[h1,om] = freqz(b,a,w);
m = 20*log10(abs(h1));
an = angle(h1);
%% Plotting
% MAGNITUDE PLOT
subplot(2,1,1);
plot(om/pi,m),xlabel('NORMALISED FREQUENCY'),ylabel('GAIN IN
DB');
title('MAGNITUDE RESPONSE');
grid on;
% PHASE PLOT
subplot(2,1,2);
plot(om/pi,an),xlabel('NORMALISED FREQUENCY'),ylabel('PHASE');
title('PHASE RESPONSE');
grid on;

69
CALCULATIONS:

70
71
72
OUTPUT:
Output of Bilinear transformation LPF:
Enter Filter Specifications
Enter Sampling frequency in Hz : 7000
Stopband frequency in Hz : 3000
Stopband attenuation in dB : 60
Passband frequency in Hz : 1500
Passband attenuation in dB : 0.15

b=
1.0e-04 *
Columns 1 through 10

0.0001 0.0014 0.0077 0.0258 0.0580 0.0929 0.1084 0.0929 0.0580 0.0258
Columns 11 through 13

0.0077 0.0014 0.0001

a=
Columns 1 through 10

1.0000 -8.1483 30.8855 -71.8971 114.3341 -130.7166 110.0744 -68.7382 31.5722 -


10.3962
Columns 11 through 13

2.3284 -0.3183 0.0201

73
74
Output of Bilinear transformation HPF:
Enter Filter Specifications
Enter Sampling frequency in Hz : 7000
Stopband frequency in Hz : 1500
Stopband attenuation in dB : 60
Passband frequency in Hz : 3000
Passband attenuation in dB : 0.15

b=
Columns 1 through 10

0.0541 -0.6491 3.5700 -11.9000 26.7751 -42.8401 49.9801 -42.8401 26.7751 -


11.9000
Columns 11 through 13

3.5700 -0.6491 0.0541

a=
Columns 1 through 10

1.0000 -6.3474 19.2890 -36.7827 48.7505 -47.1222 33.9598 -18.3418 7.3543 -


2.1314
Columns 11 through 13

0.4233 -0.0516 0.0029

75
76
Output of Impulse invariance method (LPF):
Enter Filter Specifications
Enter Sampling frequency in Hz : 7000
Stopband frequency in Hz : 3000
Stopband attenuation in dB : 60
Passband frequency in Hz : 1500
Passband attenuation in dB : 0.15

b=
1.0e-05 *
Columns 1 through 10

0.0000 -0.0001 0.0005 0.0039 0.0776 0.3522 0.6283 0.4536 0.1358 0.0151
Columns 11 through 14

0.0006 -0.0000 0.0000 0

a=
Columns 1 through 10

1.0000 -8.8797 36.8719 -94.6813 167.5781 -215.6989 207.5743 -151.1175 83.1822 -


34.1753
Columns 11 through 14

10.1821 -2.0824 0.2619 -0.0153

77
78
Output of Impulse invariance method (HPF):
Enter Filter Specifications
Enter Sampling frequency in Hz : 7000
Stopband frequency in Hz : 1500
Stopband attenuation in dB : 60
Passband frequency in Hz : 3000
Passband attenuation in dB : 0.15

b=
Columns 1 through 10

-6.0487 43.9138 -152.5072 333.4765 -509.7406 572.4557 -483.3529 308.5999 -147.5927


51.4719
Columns 11 through 14

-12.4005 1.8493 -0.1289 -0.0000

a=
Columns 1 through 10

1.0000 -7.1303 24.2875 -52.0379 77.9411 -85.9203 71.5528 -45.4913 22.0453 -


8.0318
Columns 11 through 14

2.1360 -0.3922 0.0445 -0.0024

79
80
RESULT:

81
DESIGN OF FIR DIGITAL FILTERS USING RECTANGLE WINDOW
AIM: To design and implementation of FIR filter (Low pass\High pass\Band pass\Band stop) to meet the
given specifications using Rectangular window.

SOFTWARE REQUIRED: MATLAB

THEORY:
A Finite Impulse Response (FIR) filter is a discrete linear time-invariant system whose output is based on
the weighted summation of a finite number of past inputs. A FIR transversal filter structure can be obtained
directly from the equation for discrete-time convolution.
𝑁−1

𝑦 (𝑛 ) = ∑ 𝑥 ( 𝑘 )ℎ(𝑛 − 𝑘 ) , 0 < 𝑛 < 𝑁 − 1


𝑘=0
In this equation, x(k) and y(n) represent the input to and output from the filter at time n. h(n-k) is the
transversal filter coefficients at time n. These coefficients are generated by using FDS (Filter Design Software or
Digital filter design package).
FIR – filter is a finite impulse response filter. Order of the filter should be specified. Infinite response is
truncated to get finite impulse response. placing a window of finite length does this. Types of windows available
are Rectangular, Barlett, Hamming, Hanning, Blackmann window etc. This FIR filter is an all zero filter.
The Rectangular window is given by
1, 𝑓𝑜𝑟 − (𝑁 − 1)/2 ≤ 𝑛 ≤ (𝑁 − 1)/2
𝑤𝑅 (𝑛) = {
0, 𝑂𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒

ALORITHM:
Step-1: Enter the passband ripple (𝑟𝑝 ) and stopband ripple (𝑟𝑠 ).

Step-2: Enter the passband frequency (𝑓𝑝 ) and stopband frequency (𝑓𝑠 ).

Step-3: Get the sampling frequency (𝑓).


Step-4: Calculate the normalized passband, stopband frequencies i.e. 𝑤𝑝 , 𝑤𝑠 respectively.

𝑤𝑝 = 2𝑓𝑝 ⁄𝑓 ; 𝑤𝑠 = 2𝑓𝑠 ⁄𝑓

Step-5: Calculate the order of filter using the following formula

(−20 𝑙𝑜𝑔 √𝑟𝑝 𝑟𝑠 ) − 13


𝑛 =
(14.6(𝑓𝑠 − 𝑓𝑝 )/𝑓)

[use ‘ceil()’ function for rounding off the value of n to nearest integer]
If n is an odd number then reduce its value by ‘1’.
Step-6: Generate (𝑛 + 1)𝑡ℎ point window coefficients using y = boxcar(𝑛 + 1) function for rectangle window.
Step-7: Design an 𝑛𝑡ℎ order FIR filter using the previously generated (𝑛 + 1) length window function b =
fir1(𝑛, 𝑤𝑝 , 𝑦)

82
Step-8: Find the frequency response using ‘freqz( )’ function.
Step-9: Calculate and plot the magnitude response of filter.

PROCEDURE:
1. Start the MATLAB program.
2. Create a NEW M-File from the file Menu.
3. An editor window open, type the program in it.
4. Save the file 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.

SOURCE CODE:
%% FIR filter Rectangular window
clc;
clear all;
close all;
%% Filter Specifications
disp('Enter Filter Specifications');
f = input('Enter Sampling frequency in Hz : ');
% inputs for stop band
rs = input('Stopband attenuation in dB : ');
fs = input('Stopband frequency in Hz : ');
ws = 2*fs/f;
% inputs for pass band
rp = input('Passband attenuation in dB : ');
fp = input('Passband frequency in Hz : ');
wp = 2*fp/f;
%% formulae
num = -20*(log10(sqrt(rp*rs)))-13;
dem = 14.6 * ((fs-fp)/f);
n = ceil(num/dem);
if(rem(n,2) ~= 0)
n = n-1;
end
n1 = n+1;
%% Filter Design
y = boxcar(n1);
%% LPF
bLPF = fir1(n,wp,y)
[h,o] = freqz(bLPF,1,256);
mLPF = 20*(log10(abs(h)));
%% HPF
bHPF = fir1(n,wp,'high',y)
[h,o] = freqz(bHPF,1,256);

83
mHPF = 20*(log10(abs(h)));
%% BPF
wn = [wp ws];
bBPF = fir1(n,wn,y)
[h,o] = freqz(bBPF,1,256);
mBPF = 20*(log10(abs(h)));
%% BSF
wn = [wp ws];
bBSF = fir1(n,wn,'stop',y)
[h,o] = freqz(bBSF,1,256);
mBSF = 20*(log10(abs(h)));
%% plotting
% Low Pass Filter
figure;
plot(o/pi,mLPF);
xlabel('NORMALISED FREQUENCY --->');
ylabel('Gain in dB --->');
title('Low Pass Filter');
% High Pass Filter
figure;
plot(o/pi,mHPF);
xlabel('NORMALISED FREQUENCY --->');
ylabel('Gain in dB --->');
title('High Pass Filter');
% Band Pass Filter
figure;
plot(o/pi,mBPF);
xlabel('NORMALISED FREQUENCY --->');
ylabel('Gain in dB --->');
title('Band Pass Filter');
% Band stop Filter
figure;
plot(o/pi,mBSF);
xlabel('NORMALISED FREQUENCY --->');
ylabel('Gain in dB --->');
title('Band Stop Filter');

84
OUTPUT:
Enter Filter Specifications
Enter Sampling frequency in Hz : 10000
Stopband attenuation in dB : 0.01
Stopband frequency in Hz : 1500
Passband attenuation in dB : 0.02
Passband frequency in Hz : 1000

85
bLPF =

Columns 1 through 10
-0.0112 0.0000 0.0128 0.0224 0.0242 0.0163 -0.0000 -0.0200 -0.0364 -0.0416
Columns 11 through 20
-0.0300 0.0000 0.0450 0.0970 0.1455 0.1798 0.1922 0.1798 0.1455 0.0970
Columns 21 through 30
0.0450 0.0000 -0.0300 -0.0416 -0.0364 -0.0200 -0.0000 0.0163 0.0242 0.0224
Columns 31 through 33
0.0128 0.0000 -0.0112

bHPF =

86
Columns 1 through 10
0.0115 0.0000 -0.0131 -0.0229 -0.0248 -0.0167 -0.0000 0.0205 0.0372 0.0425
Columns 11 through 20
0.0307 0.0000 -0.0460 -0.0993 -0.1489 -0.1841 0.7870 -0.1841 -0.1489 -0.0993
Columns 21 through 30
-0.0460 0.0000 0.0307 0.0425 0.0372 0.0205 -0.0000 -0.0167 -0.0248 -0.0229
Columns 31 through 33
-0.0131 0.0000 0.0115

bBPF =

Columns 1 through 10
0.0203 0.0184 -0.0000 -0.0267 -0.0437 -0.0350 0.0000 0.0428 0.0656 0.0497

87
Columns 11 through 20
-0.0000 -0.0552 -0.0811 -0.0590 -0.0000 0.0610 0.0867 0.0610 -0.0000 -0.0590
Columns 21 through 30
-0.0811 -0.0552 -0.0000 0.0497 0.0656 0.0428 0.0000 -0.0350 -0.0437 -0.0267
Columns 31 through 33
-0.0000 0.0184 0.0203

bBSF =

Columns 1 through 10
-0.0234 -0.0212 -0.0000 0.0309 0.0505 0.0404 -0.0000 -0.0494 -0.0757 -0.0573
Columns 11 through 20
0.0000 0.0637 0.0936 0.0682 0 -0.0705 0.9006 -0.0705 0 0.0682

88
Columns 21 through 30
0.0936 0.0637 0.0000 -0.0573 -0.0757 -0.0494 -0.0000 0.0404 0.0505 0.0309
Columns 31 through 33
-0.0000 -0.0212 -0.0234

RESULT:

89
DESIGN OF FIR DIGITAL FILTERS USING TRIANGLE WINDOW
AIM: To design and implementation of FIR filter (Low pass\High pass\Band pass\Band stop) to meet the
given specifications using Triangular window.

SOFTWARE REQUIRED: MATLAB

THEORY:
A Finite Impulse Response (FIR) filter is a discrete linear time-invariant system whose output is based on
the weighted summation of a finite number of past inputs. A FIR transversal filter structure can be obtained
directly from the equation for discrete-time convolution.
𝑁−1

𝑦 (𝑛 ) = ∑ 𝑥 ( 𝑘 )ℎ(𝑛 − 𝑘 ) , 0 < 𝑛 < 𝑁 − 1


𝑘=0
In this equation, x(k) and y(n) represent the input to and output from the filter at time n. h(n-k) is the
transversal filter coefficients at time n. These coefficients are generated by using FDS (Filter Design Software or
Digital filter design package).
FIR – filter is a finite impulse response filter. Order of the filter should be specified. Infinite response is
truncated to get finite impulse response. placing a window of finite length does this. Types of windows available
are Rectangular, Barlett, Hamming, Hanning, Blackmann window etc. This FIR filter is an all zero filter.does
this. Types of windows available are Rectangular, Barlett, Hamming, Hanning, Blackmann window etc. This FIR
filter is an all zero filter.
|𝑛|
The triangular window may be defined by 𝑤(𝑛) = 𝑤𝑅 [1 − 𝑀−1⁄2] ; 𝑛𝜖 [−𝑀−1 2
,
𝑀−1
2
]

ALORITHM:
Step-1: Enter the passband ripple (𝑟𝑝 ) and stopband ripple (𝑟𝑠 ).

Step-2: Enter the passband frequency (𝑓𝑝 ) and stopband frequency (𝑓𝑠 ).

Step-3: Get the sampling frequency (𝑓).


Step-4: Calculate the normalized passband, stopband frequencies i.e. 𝑤𝑝 , 𝑤𝑠 respectively.

𝑤𝑝 = 2𝑓𝑝 ⁄𝑓 ; 𝑤𝑠 = 2𝑓𝑠 ⁄𝑓

Step-5: Calculate the order of filter using the following formula

(−20 𝑙𝑜𝑔 𝑟𝑝 𝑟𝑠 ) − 13
𝑛 =
(14.6(𝑓𝑠 − 𝑓𝑝 )/𝑓)

[use ‘ceil()’ function for rounding off the value of n to nearest integer]
If n is an odd number then reduce its value by ‘1’.
Step-6: Generate (𝑛 + 1)𝑡ℎ point window coefficients using y = triang(𝑛 + 1) function for triangle window.
Step-7: Design an 𝑛𝑡ℎ order FIR filter using the previously generated (𝑛 + 1) length window function b =
fir1(𝑛, 𝑤𝑝 , 𝑦)

90
Step-8: Find the frequency response using ‘freqz( )’ function.
Step-9: Calculate and plot the magnitude response of filter.

PROCEDURE:
1. Start the MATLAB program.
2. Create a NEW M-File from the file Menu.
3. An editor window open, type the program in it.
4. Save the file 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.

SOURCE CODE:
%% FIR filter Triangle window
clc;
clear all;
close all;
%% Filter Specifications
disp('Enter Filter Specifications');
f = input('Enter Sampling frequency in Hz : ');
% inputs for stop band
rs = input('Stopband attenuation in dB : ');
fs = input('Stopband frequency in Hz : ');
ws = 2*fs/f;
% inputs for pass band
rp = input('Passband attenuation in dB : ');
fp = input('Passband frequency in Hz : ');
wp = 2*fp/f;
%% formulae
num = -20*(log10(sqrt(rp*rs)))-13;
dem = 14.6 * ((fs-fp)/f);
n = ceil(num/dem);
if(rem(n,2) ~= 0)
n = n-1;
end
n1 = n+1;
%% Filter Design
y = triang(n1);
% LPF
bLPF = fir1(n,wp,y)
[h,o] = freqz(bLPF,1,256);
mLPF = 20*(log10(abs(h)));
% HPF
bHPF = fir1(n,wp,'high',y)
[h,o] = freqz(bHPF,1,256);

91
mHPF = 20*(log10(abs(h)));
% BPF
wn = [wp ws];
bBPF = fir1(n,wn,y)
[h,o] = freqz(bBPF,1,256);
mBPF = 20*(log10(abs(h)));
% BSF
wn = [wp ws];
bBSF = fir1(n,wn,'stop',y)
[h,o] = freqz(bBSF,1,256);
mBSF = 20*(log10(abs(h)));
%% plotting
% Low Pass Filter
figure;
plot(o/pi,mLPF);
xlabel('NORMALISED FREQUENCY --->');
ylabel('Gain in dB --->');
title('Low Pass Filter');
% High Pass Filter
figure;
plot(o/pi,mHPF);
xlabel('NORMALISED FREQUENCY --->');
ylabel('Gain in dB --->');
title('High Pass Filter');
% Band Pass Filter
figure;
plot(o/pi,mBPF);
xlabel('NORMALISED FREQUENCY --->');
ylabel('Gain in dB --->');
title('Band Pass Filter');
% Band stop Filter
figure;
plot(o/pi,mBSF);
xlabel('NORMALISED FREQUENCY --->');
ylabel('Gain in dB --->');
title('Band Stop Filter');

92
OUTPUT:
Enter Filter Specifications
Enter Sampling frequency in Hz : 10000
Stopband attenuation in dB : 0.01
Stopband frequency in Hz : 1500
Passband attenuation in dB : 0.02
Passband frequency in Hz : 1000

93
bLPF =

Columns 1 through 10
-0.0007 0.0000 0.0025 0.0058 0.0078 0.0063 -0.0000 -0.0103 -0.0211 -0.0269
Columns 11 through 20
-0.0213 0.0000 0.0378 0.0877 0.1410 0.1859 0.2111 0.1859 0.1410 0.0877
Columns 21 through 30
0.0378 0.0000 -0.0213 -0.0269 -0.0211 -0.0103 -0.0000 0.0063 0.0078 0.0058
Columns 31 through 33
0.0025 0.0000 -0.0007

94
bHPF =
Columns 1 through 10
0.0007 0.0000 -0.0024 -0.0055 -0.0075 -0.0060 -0.0000 0.0098 0.0201 0.0256
Columns 11 through 20
0.0203 0.0000 -0.0360 -0.0836 -0.1343 -0.1771 0.8044 -0.1771 -0.1343 -0.0836
Columns 21 through 30
-0.0360 0.0000 0.0203 0.0256 0.0201 0.0098 -0.0000 -0.0060 -0.0075 -0.0055
Columns 31 through 33
-0.0024 0.0000 0.0007

95
bBPF =
Columns 1 through 10
0.0019 0.0035 -0.0000 -0.0103 -0.0210 -0.0202 0.0000 0.0329 0.0567 0.0477
Columns 11 through 20
-0.0000 -0.0636 -0.1012 -0.0794 -0.0000 0.0938 0.1415 0.0938 -0.0000 -0.0794
Columns 21 through 30
-0.1012 -0.0636 -0.0000 0.0477 0.0567 0.0329 0.0000 -0.0202 -0.0210 -0.0103
Columns 31 through 33
-0.0000 0.0035 0.0019

96
bBSF =
Columns 1 through 10
-0.0014 -0.0025 -0.0000 0.0074 0.0151 0.0145 -0.0000 -0.0236 -0.0407 -0.0343
Columns 11 through 20
0.0000 0.0457 0.0727 0.0570 0 -0.0674 0.9151 -0.0674 0 0.0570
Columns 21 through 30
0.0727 0.0457 0.0000 -0.0343 -0.0407 -0.0236 -0.0000 0.0145 0.0151 0.0074
Columns 31 through 33
-0.0000 -0.0025 -0.0014

97
RESULT:

98
CYCLE - 3

99
Histogram Equalization
Aim: To obtain histogram equalization of an image.
Software Required: MATLAB
Theory:
Histogram equalization is a method in image processing of contrast adjustment using
the image’s histogram. Histogram equalization often processing unrealistic effects in
photographs; however, it is very useful for scientific images like thermal or X-ray images, often
the same class of images to which one would apply false color. Also, histogram equalization
can produce undesirable effects when applied to images with low color depth. For example, if
applied to 8-bit image displayed with 8-bit grayscale it will further reduce color depth (number
of unique shades of gray) of the image. Histogram equalization will work the best when applied
to images with much higher depth than palette size, like continuous data or 16-bit grayscale
images.

Algorithm:
1. Read the image by using ‘imread()’ command.
2. Perform equalization operation on image using ‘histeq()’ command. For color image
splot the RBG planes and perform the histogram equalization on each layer.
3. Plot the images and there histograms by using ‘imshow()’ and ‘imhist()’ respectively.

Procedure:
1. Start the MATLAB program.
2. Create a NEW M-File from the file Menu.
3. An editor window open, type the program in it.
4. Save the file 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.

Source Code:
clc;
clear all;
close all;
%% Histogram equalization of grayscale image
% Read an image
a=imread('cameraman.tif');
% Histogram equalization
b=histeq(a);

100
% Plotting images and their histograms
figure;
subplot(2,2,1);imshow(a);title('Original image');
subplot(2,2,2);imshow(b);title('Histogram
equalized');
subplot(2,2,3);imhist(a);title('Histogram of original
image')
subplot(2,2,4);imhist(b);title('Histogram
equalized');
%% histogram equalization of color image
% Reading an color image
is=imread('onion.png');
% RBG plane of an image
isr=is(:,:,1); % Red plane
isg=is(:,:,2); % Green plane
isb=is(:,:,3); % Blue plane
% Plotting the images
figure,
subplot(2,2,1);imshow(is) ;title('Original image');
subplot(2,2,2);imshow(isr);title('Red plane');
subplot(2,2,3);imshow(isg);title('Green plane');
subplot(2,2,4);imshow(isb);title('Blue plane');
% Histogram equalization of an image
idr=histeq(isr);
idg=histeq(isg);
idb=histeq(isb);
% Histogram equalized image
id(:,:,1)=idr;
id(:,:,2)=idg;
id(:,:,3)=idb;
% Plotting the images
figure,
subplot(2,2,1);imshow(id);title('Histogram equalized
image');

101
subplot(2,2,2);imshow(idr);title('Red plane');
subplot(2,2,3);imshow(idg);title('Green plane');
subplot(2,2,4);imshow(idb);title('Blue plane');
figure,
subplot(2,2,1);imshow(is);title('Original image');
subplot(2,2,2);imhist(isr);title('Histogram of Red');
subplot(2,2,3);imhist(isg);title('Histogram of
Green');
subplot(2,2,4);imhist(isb);title('Histogram of
Blue');
figure,
subplot(2,2,1);imshow(id);title('Histogram equalized
image');
subplot(2,2,2);imhist(idr), title('Histogram of
Red');
subplot(2,2,3);imhist(idg),title('Histogram of
Green');
subplot(2,2,4);imhist(idb), title('Histogram of
Blue');

102
Output:
Grayscale image –

a. Histogram equalization of given grayscale image

103
104
Colour image –

i. RGB planes of a given colour image

ii. RGB plains of a Histogram equalized image

105
iii. Histogram of RGB plains of given image

iv. Histogram of RGB plains of Histogram equalized image

106
Result:

107
Edge Detection
Aim: To detect the edge of the grayscale images using with and without command.
Software Required: MATLAB
Theory:
Edges characterize boundaries and are therefore a problem of fundamental importance
in image processing. Edges in in=mages are areas with strong intensity contrast – a jump in
intensity from one pixel to the next. Edge detecting an image significantly reduces the amount
of data and filters out useless information, while preserving the important structural properties
in an image. There are many ways to perform edge detection. However, the majority of
different methods may be grouped into two categories, Gradient and Laplacian. The Gradient
method detects the edges by looking for the maximum and minimum in the first derivative of
the image. The Laplacian method searches for zero crossing in the second derivative of the
image to find edges. An edge has the one-dimensional shape of a ramp and calculating the
derivative of the image can highlight its location.

Algorithm:
1. Read the image by using ‘imread()’ command.
2. Perform edge operation using ‘edge(I, method)’ where I is given image; methods are
‘sobel’, ‘roberts’ and ‘perwett’.
3. Plot the image using ‘imshow()’ command.

Procedure:
1. Start the MATLAB program.
2. Create a NEW M-File from the file Menu.
3. An editor window open, type the program in it.
4. Save the file 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.

Source Code:
clc;
clear all;
close all;
%% First order derivatives
a=imread('cameraman.tif');
e1=edge(a,'sobel');
e2=edge(a,'roberts');
e3=edge(a,'prewitt');
figure,

108
subplot(2,2,1);imshow(a); title('Original image');
subplot(2,2,2);imshow(e1);title('Sobel');
subplot(2,2,3);imshow(e2);title('Roberts');
subplot(2,2,4);imshow(e3);title('Prewitt');
%% without edge command
i=imread('cameraman.tif');
p=[0,1,0; 1,-4,1;0,1,0];
icx=filter2(p,i);
py=p';
icy=filter2(py,i);
pedge=sqrt(icx.^2 + icy.^2);
figure,
subplot(2,2,1);imshow(i);title('Original image')
subplot(2,2,2);imshow(uint8(pedge));title('Laplacian')
%arithmetic operations
a=uint8(i)-uint8(pedge);
subplot(2,2,3);imshow(a);title('Image subtraction')
b=uint8(i)+uint8(pedge);
subplot(2,2,4);imshow(b);title('Image addition')

109
Output:
With using Function

110
Without using fuction

111
Result:

112

You might also like