DSP Lab Manual final
DSP Lab Manual final
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.
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
Exponential sequence:
𝑥 (𝑛 ) = 𝑎 𝑛
Sinusoidal sequence:
𝑥(𝑛) = 𝐴 sin(𝜔𝑛 + ф)
Square sequence:
𝑁
1, 0≤𝑛≤
𝑥 (𝑛 ) = { 2
𝑁
−1, <𝑛≤𝑁
2
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
8
Unit Ramp function:
Exponential function:
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.
THEORY:
The mathematical definition of convolution in discrete time domain is
∞
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)');
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.
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
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=
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.
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
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
40
enter modulation index: 0.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.
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 (𝑓𝑠 ).
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’)
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;
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.
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 (𝑓𝑠 ).
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;
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;
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;
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=
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=
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=
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.
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 (𝑓𝑠 ).
64
Step-3: Get the sampling frequency (𝑓).
Step-4: Calculate the normalized passband, stopband frequencies i.e. 𝑤𝑝 , 𝑤𝑠 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:
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;
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;
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;
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
a=
Columns 1 through 10
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
a=
Columns 1 through 10
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
a=
Columns 1 through 10
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
a=
Columns 1 through 10
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.
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
ALORITHM:
Step-1: Enter the passband ripple (𝑟𝑝 ) and stopband ripple (𝑟𝑠 ).
Step-2: Enter the passband frequency (𝑓𝑝 ) and stopband frequency (𝑓𝑠 ).
𝑤𝑝 = 2𝑓𝑝 ⁄𝑓 ; 𝑤𝑠 = 2𝑓𝑠 ⁄𝑓
[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.
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
ALORITHM:
Step-1: Enter the passband ripple (𝑟𝑝 ) and stopband ripple (𝑟𝑠 ).
Step-2: Enter the passband frequency (𝑓𝑝 ) and stopband frequency (𝑓𝑠 ).
𝑤𝑝 = 2𝑓𝑝 ⁄𝑓 ; 𝑤𝑠 = 2𝑓𝑠 ⁄𝑓
(−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 –
103
104
Colour image –
105
iii. Histogram of RGB plains of given 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