Digital signal processing laboratory (18EEL67)

STUDY EXPERIMENT
1. Program to generate a triangular and Square wave using Matlab.
%Generation of Square wave
A2=5;
f2=10;
t2=linspace(0,1,1001);
y2=A2*square(2*pi*f2*t2,50);
fh2=figure('Name','Square wave of 50% duty cycle');
ah2=axes;
ph2=plot(t2,y2)
xlabel(ah2,'Time in seconds');
ylabel(ah2,'Amplitude');
title(ah2,'Amplitude of square wave vs time');
set(ah2,'Ylim',[-10 10])

%Generation of Triangular Wave

A3=12;%
f3=12;
t3=linspace(0,1,1001);
y3=A3*sawtooth(2*pi*f3*t3,0.5);
fh3=figure('Name','Triangular Wave');
ah3=axes;
ph3=plot(t3,y3)
set(ah3,'Ylim',[-12 12])
xlabel(ah3,'Time in seconds');
ylabel(ah3,'Amplitude');
title(ah3,'Amplitude of Triangular wave vs time');
Amplitude of square wave vs time
10

2
Amplitude

-2

-4

-6

-8

-10
0 0.2 0.4 0.6 0.8 1
Time in seconds

Amplitude of Triangular wave vs time

10

5
Amplitude

-5

-10

0 0.2 0.4 0.6 0.8 1

Time in seconds

Dept. of EEE, Citech Page 1

 Digital signal processing laboratory (18EEL67)

2. Program to generate a unit circle of given radius.

Theta=linspace(0,2*pi,100);
X=cos(theta);
Y= sin(theta);
Plot(x,y)
Axis(‘equal’);
Xlabel(‘x’)
Ylabel(‘y’)
Title(‘circle of unit radius’)
print

Dept. of EEE, Citech Page 2

 Digital signal processing laboratory (18EEL67)

EXPERIMENT No 1: Verification of Sampling Theorem.

Aim : To verify Sampling theorem for a signal of given frequency.

Theory:
• Sampling is a process of converting a continuous time signal (analog signal) x(t) into a discrete
time signal x[n], which is represented as a sequence of numbers. (A/D converter)

• Converting back x[n] into analog (resulting in x (t ) ) is the process of reconstruction. (D/A
converter)
• Sampling Theorem can be stated in two forms-
o If a finite energy signal g(t) contains no frequencies higher than ‘w’ Hz, it is completely
determined by specifying its ordinates at a sequence of points spaced 1/2w sec apart.
o If a finite energy signal g(t) contains no frequencies higher than ‘w’ Hz, it may be
completely recovered from its ordinates at a sequence of points spaced 1/2w sec apart.
• The sampling frequency f s determine the spacing between samples. The minimum sampling
rate of 2w samples per second, for a signal bandwidth of ‘w’ Hz is called NYQUIST RATE.
ƒs = 2ƒm where,
ƒs=Sampling Frequency, ƒm=Message frequency
We encounter these conditions in Sampling Theorem, they are Right sampling, Under sampling
and Over sampling
➢ Nyquist Sampling → ƒs = 2ƒm
➢ Under Sampling → ƒs ≤ 2ƒm
➢ Over Sampling → ƒs ≥ 2ƒm
• Under sampling:
When the sampling rate is less than the nyquist rate, there is overlap of the frequency
components and the higher frequency components and the other higher frequency terms of the
continuous time signals appear in the low frequency region. This is referred to as aliasing.
• Over sampling:
If the sampling rate is greater then there is a lot of redundancy in the discrete time signal. This
is referred to as over sampling.

MATLAB Implementation:
Step 1: a) Generate a sinusoidal signal of length 0. 05secs.For an approximate analog signal xt
,choose the spacing between the samples to be very small (≈0), say 50µs = 0.00005 [Smaller the
spacing, better will be the approximation to analog signal].
b) The vector that represents the time base t = 0:0.00005:0.05;
The colon (:) operator in MATLAB creates a vector, in the above case a time vector running from
0 to 0.05 in steps of 0.00005 i.e., 1000 sample points are created. The semicolon (;) tells MATLAB
not to display the result.
c) Given t, the analog signal xt of desired frequency fd is generated as (sin(ωt) = sin(2πft))
xt=sin(2*pi*fd*t);
pi is recognized as 3.14 by MATLAB.
Step 2: To illustrate under sampling, nyquist sampling and oversampling condition, choose
sampling frequency fs1<2fd, fs2=2fd and fs3>>2fd. For this sampling rate T*=1/fs*

Algorithm:

Dept. of EEE, Citech Page 3

 Digital signal processing laboratory (18EEL67)

1. Input the desired frequency fd (for which sampling theorem is to be verified).

2. Generate an analog signal xt of frequency fd.
3. Generate under sampled, nyquist and oversampled discrete time signals.
4. Plot the waveforms and verify the sampling theorem.

MATLAB Program:

tfinal=0.05;
t=0:0.00005:tfinal; % define time vector .05/.00005=1000 sample pts
fd=input('Enter analog frequency'); %200 Hz
xt=sin(2*pi*fd*t); %Generate analog signal
%condition for under sampling i.e., fs1<2*fd
fs1=1.3*fd;
n1=0:1/fs1:tfinal; % define time vector .05/.00384=13 sample pts
xn=sin(2*pi*n1*fd); %Generate under sampled signal
subplot(3,1,1);
plot(t,xt,'b',n1,xn,'r*-');
title('under sampling plot');

%condition for Nyquist sampling fs2=2*fd

fs2=2*fd;
n2=0:1/fs2:tfinal; %define the time vector 05/.00384=13 sample pts
xn=sin(2*pi*fd*n2); % Generate sampled signal
subplot(3,1,2);
plot(t,xt,'b',n2,xn,'r*-');
title('Nyquist plot');

%condition for oversampling i.e., fs1>2*fd

fs3=5*fd;
n3=0:1/fs3:tfinal; %define the time vector 05/.001=50 sample pts
xn=sin(2*pi*fd*n3); % Generate over sampled signal
subplot(3,1,3);
plot(t,xt,'b',n3,xn,'r*-');
title('Oversampling plot');
xlabel('time');
ylabel('amplitude');
legend('analog','discrete')

Inference:
1. From the under-sampling plot observe the aliasing effect. The analog signal is of 200Hz
(T=0.005s). The reconstructed (from under sampled plot) is of a lower frequency. The alias
frequency is computed as
fd-fs1 = 200-1.3*200 = 200-260= -60Hz
This is verified from the plot. The minus sign results in a 180˚ phase shift.
(For example: 3kHz & 6kHz sampled at 5kHz result in aliases of -2kHz (3k-5k) & 1kHz
(6k-5k) respectively)

Dept. of EEE, Citech Page 4

 Digital signal processing laboratory (18EEL67)

2. Sampling at the Nyquist rate results in samples sin(πn) which are identically zero, i.e., we
are sampling at the zero crossing points and hence the signal component is completely
missed. This can be avoided by adding a small phase shift to the sinusoid. The above
problem is not seen in cosine waveforms (except cos(90n)). A simple remedy is to sample
the analog signal at a rate higher than the Nyquist rate. The Fig1.2 shows the result due to
a cosine signal (x1=cos(2*pi*fd*n1);
3. The over sampled plot shows a reconstructed signal almost similar to that of the analog
signal. Using low pass filtering the wave form can be further smoothened.

Fig.1.1 Plots of a sampled sine wave of 200Hz

Fig.1.2 Plots of a sampled cosine wave of 200Hz

Assignment:
1. Generate sine/cosine waveforms of multiple frequencies, say x1=sin(2*pi*fd1*n1)+
sin(2*pi*fd2*n2); where fd1, fd2 are 2 different frequencies. Verify sampling theorem.

Dept. of EEE, Citech Page 5

 Digital signal processing laboratory (18EEL67)

Viva Questions:
• State sampling theorem?
• What do you mean by process of reconstruction.
• What do you mean Aliasing? What is the condition to avoid aliasing for sampling?
• Write the conditions of sampling.
• Explain the statement t= 0:0.000005:0.05
a. What does colon (:) and semicolon (;) denotes.
• What is the relation between the continuous time angular frequency and normalized discrete
time angular frequency?
• What is the relation between the bandwidth of continuous time signal and the bandwidth of
corresponding discrete time signal after sampling process?
• What is Nyquist frequency, Nyquist rate, Nyquist band?
• What is the bandwidth of the continuous time signal?
• What is the use of command 'legend'?

EXPERIMENT NO. 2: IMPULSE RESPONSE OF A GIVEN SYSTEM

Dept. of EEE, Citech Page 6

 Digital signal processing laboratory (18EEL67)

Aim: To find the impulse response h(n) of the LTI system whose response y(n) to an input x(n) is
given.
Theory:

Fig.2.1 A LTI system

• A discrete time LTI system (also called digital filters) is represented by

o A linear constant coefficient difference equation
a0y[n] + a1y[n - 1] + …+ amy[n - m] = b0x[n] + b1x[n - 1] + b2x[n - 2] + …+ bkx[n - k]
o System function H(z) (obtained by applying Z transform to the difference equation).
Y ( z ) b0 + b1 z −1 + b2 z −2 + ....bk z − k
H ( z) = =
X ( z ) 1 + a1 z −1 + a2 z −2 + ....am z −m
The transfer function is defined in terms of the coefficients of the numerator polynomial
and the denominator polynomial. The inverse Z-transform of the transfer function gives
the impulse response of the system.

o A frequency response function obtained by applying DTFT on the impulse response

− j −2 j
h[n] (or by replacing z with ejΩ in H(z)) to get H (e j ) = Y (e ) = b0 + b1e + b2 e
j

j − j − 2 j
X (e ) 1 + a1e + a2 e

MATLAB Implementation:
• Given the input sequence x[n] and the output sequence y[n], we can find the impulse function
h[n] by using the inverse operation deconv (program 1).
h=deconv(y,x)
Note: y the output sequence should be of longer length than x.
➢ The deconvolution operation is valid only if the LTI system is ‘invertible’.
➢ Valid y and x should be given, else the conv verification will result in a slightly
different y (because of the remainder generated by deconv).
• If both y[n] and x[n] are finite then the impulse response is finite (FIR).
• Given the difference equation or H(z), the impulse response of the LTI system is found using
filter or impz MATLAB functions (program 2).
Problem Statement 1
Find the impulse response h(n) of the given LTI system whose response is y(n)=[1 2 1 1 0 -1]
to an input x(n)=[1 1 1 1]. Also verify the result using conv operation (FIR).

Algorithm:
1. Input the two sequences x and y.
2. Use deconv to get impulse response h.
3. Plot the sequences.
4. Verify values of y using conv(x,h).

MATLAB Programs:

Dept. of EEE, Citech Page 7

 Digital signal processing laboratory (18EEL67)

Program 1:
y= input('The output sequence y(n) of the system=');
x=input('the input sequence of the system=');
[h,r]=deconv(y,x);
disp('the impulse response of the system is=');
disp(h);
N=length(h);
n=0:1:N-1;
stem(n,h); %graphical display part
xlabel('Time index n');
ylabel('Amplitude');
title('impulse response of a system')

yv=conv(x,h)+r; %Verification
disp(‘the verified output sequence is’);
disp(yv)

Result:
The output sequence y(n) of the system=[1 2 1 1 0 -1]
the input sequence of the system=[1 1 1 1]
the impulse response of the system is= 1 1 -1
the verified output sequence is 1 2 1 1 0 -1

Inference:
The impulse response h[n] is of finite duration.
The verified convolution output sequence is the same as the given y[n].
Plot of h given x and y

Fig 2.1 Impulse Response

Problem Statement 2
A LTI system is described by the difference equation y(n) - y(n-1) + 0.9y(n-2) = x(n);
Find its impulse response h(n) at n=0,…..,100. (IIR)

Explanation:
As per the given difference equation the co-efficient of a0=1, a1=-1, a2=0.9 and b0=1

Algorithm:
1. Input the two coefficient arrays as b and a.
2. Input the length of the impulse response required as N=100.
3. Use impz function to get impulse response h.
4. Plot the impulse response roughly.

Dept. of EEE, Citech Page 8

 Digital signal processing laboratory (18EEL67)

5. Input the length of the impulse response required as N=10.

6. Use impz function to get impulse response h.
7. Plot the impulse response according to scale.

MATLAB Program:
%To find Impulse Response given the difference equation
N=input('Length of impulse response required=');
b=[1]; %x[n] coefficient
a=[1,-1,0.9]; %y[n] coefficients
h=impz(b,a,N); %impulse response
%Plot the Impulse Response Sequence
n=0:N-1; %time vector for plotting
stem(n,h);
xlabel('n');
ylabel('h(n)');
title('impulse response');

Results:
Length of impulse response required=100
*enter the x[n] coefficient array=[1]
*enter the y[n] coefficient array=[1 -1 0.9]

Inference:
The impulse response h[n] is of infinite duration. h[n] is stable (as h[n] is absolutely summable
+

 h[n]   , decaying exponential).

n = −

Plot of h given coefficients of x and y for N=100

Note:
• deconv -Deconvolution and polynomial division. [Q,R] = deconv (B,A) deconvolves vector
A out of vector B. The result is returned in vector Q and the remainder in vector R such that
B = conv(A,Q) + R. If A and B are vectors of polynomial coefficients, deconvolution is
equivalent to polynomial division. The result of dividing B by A is quotient Q and remainder
R. (Remember conv is polynomial multiplication).
• h=impz(b,a,N), returns the impulse response h[n], where b and a are coefficient arrays obtained
from difference equation, N = length of the response required.

Dept. of EEE, Citech Page 9

 Digital signal processing laboratory (18EEL67)

• IMPZ - Impulse response of digital filter. [H,T] = IMPZ(B,A) computes the impulse response
of the filter B/A choosing the number of samples for you, and returns the response in column
vector H and a vector of times (or sample intervals) in T (T = [0 1 2 ...]').
• [H,T] = IMPZ(B,A,N) computes N samples of the impulse response.

Do it;
1) y(n) – 1/6y(n-1) -1/6y(n-2) = x(n);
a=[1 -1/6 -1/6];
b= [1];

2) y(n) – 3/4y(n-1) +1/8y(n-2) = x(n)-1/3x(n-1);

a=[1 -3/4 1/8];
b= [1 -1/3];

Solved examples:
Example1:
Let’s take a second order difference equation
y(n) – (1/6) y(n-1) –(1/6) y(n-2) = x(n) ………………..(1)
For impulse response the particular solution yp(n) = 0 So, y(n) = yh(n) …………….(2)
Let yh(n) = λn ……………..(3)
Substituting (3) in (1) we get,
λn - (1/6) λn-1 – (1/6) λn-2 = 0 ( x(n) = 0 for homogeneous solution)
or
λ2 - (1/6) λ – (1/6) = 0
The roots of the characteristic equation are,
λ1 = 1/2, λ2 = -1/3
The solution yh(n) = C1 (1/2)n + C2 (-1/3)n………… …(4)
For impulse x(n) = δ(n); x(n) = 0 for n > 0 and x(0) = 1
Substituting the above relations in eqn (1) we have,
For n = 0,
y(0) – (1/6) y(-1) –(1/6) y(-2) = x(0) = 1
i.e. y(0) = 1. ………………….(5)
For n = 1,
y(1) – (1/6) y(0) –(1/6) y(-1) = x(1)
y(1) – (1/6) = 0
y(1) = 1/6. ……………..(6)
From eqn (4)
y(0) = C1 + C2 …………………..(7)
y(1) = (1/2)C1 – (1/3)C2 ……………(8)
Comparing equations (5), (6), (7) and (8) we get,
C1 + C2 = 1
(1/2)C1 – (1/3)C2 = (1/6)
Solving for C1 and C2 we get
C1 = 3/5, C2 = 2/5
Therefore, the solution
y(n) = (3/5)(1/2)n + (2/5)(-1/3)n …………….…(9)
Now, let’s find out impulse response.
We know that,
y(n) = x(n) * h(n)
h(n) = y(n) / x(n) …………………………..…(10)
Compute y(n) for various values of n
Substitute n = 0, 1, 2, 3 in eqn (9) we get,

Dept. of EEE, Citech Page 10

 Digital signal processing laboratory (18EEL67)

y(0) = 1.0000
y(1) = 0.1667
y(2) = 0.1944
y(3) = 0.0602
So, y(n) = {1.0000, 0.1667, 0.1944, 0.0602}
We know that, x(n) = 1
So, the impulse response for the given difference equation using eqn (10) is
h(n) = {1.0000, 0.1667, 0.1944, 0.0602}

Also slove: y(n) – (3/4) y(n-1) –(1/8) y(n-2) = x(n)-(1/3)x(n-1)

Example2: y(n)=[1 2 3 2 1 2 ] x(n) =[1 2 3] h(n)= 1 0 0 2

h(n)= [1 0 0 2] y(n) = [1 2 3 2 1 2] verified

Example 3 : y(n) =[1 2 1 1 0 -1]x(n)=[1 1 1 1] h(n)= 1 1 -1

h(n)= [1 1-1]
Viva Questions
1) What is impulse?
2) Define a unit impulse function?
3) What is meant by impulse response?
4) How can an impulse be modelled mathematically?
5) The Laplace transform of the impulse response function is known as the transfer function.
7) How sampling is done by an impulse function?
8) What is a linear system?
9) What is a time invariant system?

Dept. of EEE, Citech Page 11

 Digital signal processing laboratory (18EEL67)

EXPERIMENTNO 3: LINEAR CONVOLUTION OF TWO SEQUENCES

Aim: To obtain convolution of two finite duration sequences.

Theory:
Impulse function:
The continuous time version of the unit impulse, commonly denoted by δ(t), is defined by the
following pair of relations
δ(t) = 0 for t ≠0and
δ(t) = 1 for t = 0

The discrete-time version of the impulse response commonly denoted by δ[n], is defined by
δ[n] = 1 for n = 0;
0 for n ≠ 0

Impulse response of a linear time-invariant system:

Consider the product signal, x[n]δ[n] = x[0]δ[n]
We can generalize this relationship as x[n]δ[n-k] = x[k]δ[n-k], i.e.
we can express x[n] as a weighted sum of impulse functions as shown in the equation given
below
x[n] = ………… + x[-2] δ[n+2] + x[-1] δ[n+1] + x[0] δ[n] + x[1] δ[n-1] +….

Let H denote the system to which the input is applied. Then using the above equation to
represent the input x[n] to the system results in the output.
y[n]=H{∑x[k] δ[n-k]}
As the system is linear and time invariant, The output y[n] of a LTI (linear time invariant) system
can be obtained by convolving the input x[n] with the system’s impulse response h[n].

y[n] = ∑x[k] H{ δ[n-k]}

+ +
The convolution sum is y[n] = x[n]  h[n] =  x[k ]h[n − k ] =
k = −
 x[n − k ]h[k ]
k = −
The length of y will be length(x)+length(h)-1

MATLAB Implementation:
• If both the sequences x1 and x2 are of finite duration, then we can use the MATLAB
function ‘conv’ to evaluate the convolution & sum to obtain the output y[n]. Convolution
is implemented as polynomial multiplication (refer MATLAB help).
• The length of y[n] = xlength + hlength -1.
• The conv function assumes that the two sequences begin at n=0 and is invoked by
y=conv(x,h).
But MATLAB recognizes index 1 to positive maximum. Index 0 is not recognized.
• For two sided sequence, i.e., for sequence that exist from –n to n, say n=-3,or n=2; then it
is required to provide beginning and end point to y[n].
➢ This can be done by finding the start index of x and h and end index of x and h
and summing them up to find the start and end index of y.
➢ Note that during plotting, the time vector size and the sequence size should be the same,
i.e., the number of elements in the sequence x and in its time vector n should be the same.

Problem Statement 1
Find the convolution of two sequences given first sequence X1= [1 2 3 4] and second sequence
X2=[ 2 3 2 4]

Dept. of EEE, Citech Page 12

 Digital signal processing laboratory (18EEL67)

Algorithm:
1. Input the two sequences x1, x2
2. Convolve x1 and x2 using conv function to get output y.
3. Plot the input and output sequences.
4. Label the axis and give appropriate title.

MATLAB Program:
%Level 1 Program (both sequences are right sided)

x1=input('Enter first sequence')

x2=input('Enter second sequence')
y=conv(x1,x2);
disp('linear con of x1 & x2 is y=');
disp(y);
subplot (2,1,1);
stem(y); %graphical display part
xlabel('time index n');
ylabel('amplitude ');
title('convolution output');

subplot(2,2,3);
stem(x1);
xlabel('time index n');
ylabel('amplitude ');
title('plot of x1');

subplot(2,2,4);
stem(x2);
xlabel('time index n');
ylabel('amplitude ');
title('plot of x2');

Result
Enter first sequence [1 2 3 4]
x1 = 1 2 3 4
Enter second sequence[ 2 3 2 4]
x2 = 2 3 2 4
linear con of x1 & x2 is y=
2 7 14 25 26 20 16

verify the plots as shown in Figure3.1 window

Dept. of EEE, Citech Page 13

 Digital signal processing laboratory (18EEL67)

plot of x1 plot of x2
4 4

3 3

amplitude

amplitude
2 2

1 1

0 0
1 2 3 4 1 2 3 4
time index n time index n
convolution output
30

20
amplitude

10

0
1 2 3 4 5 6 7
time index n

Figure3.1 Linear convolution of right sided sequence

Problem Statement 2
Find the convolution of two sequences given first sequence x1=[1 2 3 2 1 3 4] & second
sequence x2=[ 2 -3 4 -1 0 1]

Algorithm:
1. Input the two sequences x1, x2 along with their time index
2. Plot the input and output sequences against their time index.
3. Convolve x1 and x2 using conv function to get output y.
4. Find the start and end time indices and define the time index for y
5. Label the axis and give appropriate title

%Level 2 Program (Two sided sequences)

x1=[1 2 3 2 1 3 4] %first sequence
n1=-3:3 %time vector of first seq
x2=[ 2 -3 4 -1 0 1] %second seq
n2=-1:4 %time vector of second seq

ybegin=n1(1)+n2(1); %add first elements of time vector

yend=n1(length(x1))+n2(length(x2)); %add last elements of time vector
ny=[ybegin:yend];

y=conv(x1,x2);
disp('linear con of x1 & x2 is y=');
disp(y);

subplot (2,1,1); stem(ny,y); %graphical display with time info

xlabel('time index n');
ylabel('amplitude ');
title('convolution output');

subplot(2,2,3);
stem(n1,x1);

Dept. of EEE, Citech Page 14

 Digital signal processing laboratory (18EEL67)

xlabel('time index n');

ylabel('amplitude ');
title('plot of x1');

subplot(2,2,4);
stem(n2,x2);
xlabel('time index n');
ylabel('amplitude ');
title('plot of x2');

Result
x1 = 1 2 3 2 1 3 4
n1 = -3 -2 -1 0 1 2 3
x2 = 2 -3 4 -1 0 1
n2 = -1 0 1 2 3 4
linear con of x1 & x2 is y=
2 1 4 2 6 9 3 2 15 -3 3 4

verify the plots as shown in Figure3.2 window

convolution output
15

10
amplitude

-5
-4 -2 0 2 4 6 8
time index n
plot of x1 plot of x2
4 4

3 2
amplitude

amplitude

2 0

1 -2

0 -4
-4 -2 0 2 4 -2 0 2 4
time index n time index n

Figure3.2 Linear convolution of both sided sequence

Dept. of EEE, Citech Page 15

 Digital signal processing laboratory (15EEL68)

Solved examples:

Example:
x(n) = {1, 2, 3, 1}
h(n) = {1, 1, 1}
length(y(n)) = length(x(n)) + length(y(n)) – 1
=4+3–1=6

Dept. of EEE, Citech Page 16

 Digital signal processing laboratory (15EEL68)

Output:
Enter the 1st seq: [1 2 3 1]
Enter the 2nd seq: [1 1 1]
The linear convolution of two sequences:
1 3 6 6 4 1

Method2: single sided sequence

y(n)= {2, 7 , 14, 25,26,20,16}

b) both sided sequence

x1 =1 2 3 2 1 3 4

x2 =2 -3 4 -1 0 1

linear convolution of x1 & x2 is y=

y(n)= 2 1 4 2 6 9 3 2 15 -3 3 4 for
n= -4 -3 -2 -1 0 1 2 3 4 5 6 7

Viva Questions:

1. Define signal and system.

2. What is linear convolution?
3. Explain how convolution syntax built in function works.
4. How to calculate the beginning and end of the sequence for the two sided controlled
output?
5. What is the total output length of linear convolution sum.
6. What is an LTI system? Write the expressions for LTI system convolution formula
7. Describe impulse response of a function.
8. Is it possible to represent any discrete time signal in terms of impulses? If yes, represent
by using example.
9. Explain scaling and superposition properties of a system.

Dept. of EEE, Citech Page 17

 Digital signal processing laboratory (15EEL68)

10. What is meant by linearity of a system and how it is related to scaling and superposition?
11. What is energy signal? How do you calculate energy of a signal?
12. What is power signal? How to calculate power of a signal?
13. Differentiate between even and odd signals.
14. Explain time invariance property of a system with an example.
15. What is memory less system?
16. When a system is said to have memory?
17. What is meant by causality?
18. Explain linear convolution and circular convolution.
19. What us the length of linear convolution if length of input & impulse responses are N1 &
N2 respectively?
20. How do you calculate the beginning and end of the sequence for the two sided controlled
output?
21. What is the length of linear and circular convolutions if the two sequences are having the
length n1 and n2?

Dept. of EEE, Citech Page 18

 Digital signal processing laboratory (15EEL68)

EXPERIMENT NO 4: Circular convolution of two given sequences.

To perform circular convolution of given sequences using (a) The convolution summation formula (b) The
matrix method and (c) Linear convolution from circular convolution with zero padding.

4a) The circular convolution using summation formula

Aim: To obtain circular convolution of two finite duration sequences using the
convolution summation formula.

Theory:
The output y[n] of a LTI (linear time invariant) system can be obtained by convolving the input
x[n] with the system’s impulse response h[n].
• Linear convolution is generally applied to aperiodic sequences. Whereas the Circular
Convolution is used to study the interaction of two signals that are periodic.
+ +
• The linear convolution sum is y[n] = x[n]  h[n] =  x[k ]h[n − k ] =
k = −
 x[n − k ]h[k ] .
k = −

• To compute circular convolution, the following operations are involved:

o Folding: fold h[k] to get h[-k] ( for y[0])
o Multiplication : vk[n] = x[k] × h[n-k] (both sequences are multiplied sample by
sample)
o Addition: Sum all the samples of the sequence vk[n] to obtain y[n]
o Shifting : the sequence h[-k] to get h[n-k] for the next n.
+
• The circular convolution sum is y[n] = x[n]( N )h[n] =  x[k ]h[ n − k N
] where the
k = −

index n − k N
implies circular shifting operation and −k N
implies folding the
sequence circularly.
• Steps for circular convolution are the same as the usual convolution, except all index
calculations are done "mod N" = "on the wheel".
o Plot f[m] and h[−m] as shown in Fig. 4.1. (used f(m) instead of x(k))
o Multiply the two sequences
o Add to get y[m]
o "Spin" h[−m] n times Anti Clock Wise (counter-clockwise) to get h[n-m].
• x[n] and h[n] can be both finite or infinite duration sequences. If infinite sequences,
they should be periodic, and the N is chosen to be at least equal to the period. If they
are finite sequences, then N is chosen as max(xlength, hlength).
• Say x[n] = {−2,3, 1,1} and N = 5, then the x[-1] and x[-2] samples are plotted at x[N-1]


= x[-4] and x[N-2] =x[-3] places. Now x[n] is entered as x[n] = {1,1,0,−2,3} .


Fig. 4.1 Plotting of f(m) and h(-m) for circular convolution

Dept. of EEE, Citech Page 19

 Digital signal processing laboratory (15EEL68)

MATLAB Implementation:
MATLAB recognizes index 1 to be positive maximum. Index 0 is not recognized.
Hence in the below program wherever y, x and h sequences are accessed, the index is added
with +1. the modulo index calculation for circular convolution is carried out using the function
MOD Modulus (signed remainder after division). MOD(x,y) is x - y.*floor(x./y) if y ~= 0.
By convention, MOD(x,0) is x. The input x and y must be real arrays of the same size, or real
scalars. MOD(x,y) has the same sign as y while REM(x,y) has the same sign as x. MOD(x,y)
and REM(x,y) are equal if x and y have the same sign, but differ by y if x and y have different
signs.

Algorithm:
1. Input the two sequences as x and h.
2. Circularly convolve both to get output y.
3. Plot the input and output sequences.

MATLAB Program:
clc
x=input('Enter the sequence 1:'); %input sequence
h=input('Enter the sequence 2:'); %impulse response
N=max(length(x),length(h)); %calculate value of N

for n=0:N-1
y(n+1)=0;
for k=0:N-1
i=mod((n-k),N); %calculation of x index
if i<0
i=i+N;
end %end of ‘if’
y(n+1)=y(n+1)+h(k+1)*x(i+1);
end %end of inner ‘for loop’
end %end of outer ‘for loop’
disp('circular convolution of x & h is y=');
disp(y);
n1=0:N-1;
stem(n1,y);
title('Circular convolution output y(n)');

Result:
Enter the sequence 1: [1 2 3 4]
Enter the sequence 2: [1 2 3 4]
circular convolution of x1 &x2 is y=
26 28 26 20

Dept. of EEE, Citech Page 20

 Digital signal processing laboratory (15EEL68)

Plot of circular convolution of x and h

4b) The circular convolution using Matrix method

Aim: To obtain circular convolution of two finite duration sequences using the matrix
method.

Theory:
The output y[n] has n varying from 0 to N-1.
N −1
The circular convolution sum is y[n] = x[n]  h[n] =  x[k ]h[ n − k N
]
k =0
Expanding for k and n from 0 to N-1, say for the specific case of N=4; the above equation in
matrix form is given below
N −1
y[n] = x[n]  h[n] =  x[k ]h[ n − k N
]
k =0
3
y[n] =  x[k ]h[ n − k 4 ]
k =0

 y[0] h[0] h[3] h[2] h[1]   x[0]

 y[1]   h[1] h[0] h[3] h[2]  x[1] 
 = 
 y[2] h[2] h[1] h[0] h[3]  x[2]
     
 y[3]  h[3] h[2] h[1] h[0]  x[3]

MATLAB:
clc;
x=input('enter the 1st sequence x(n)=')
h=input('enter the 1st sequence h(n)=')
N=max(N1,N2)
%zero padding the sequences x & h
x1=[x zeros(1,N-length(x))]
h1=[h zeros(1,N-length(h))]

x=transpose(x1)
h=transpose(h1)
temp=h;
for i=1:N-1;
temp=circshift(temp,1);
h=horzcat(h,temp);
end;
y=h*x
Dept. of EEE, Citech Page 21
 Digital signal processing laboratory (15EEL68)

disp('circular convolution of x & h is y=');

disp(y);
n=0:N-1;
stem(n,y);
title('Circular convolution output y(n)');

Result:
Enter the sequence 1: [1 2 3 4]
Enter the sequence 2: [1 2 3 4]
circular convolution of x1 &x2 is y=
26 28 26 20

Plot of circular convolution of x and h

4c) Linear convolution from circular convolution with zero padding

Aim: To obtain Linear convolution from circular convolution with zero padding

Theory:
Circular convolution is linear convolution with aliasing. Increasing the size N from
Max(xlength,hlength) to N = xlength + helngth -1 and then running the circular convolution
program results in a linear convolution output

MATLAB Program
x=input('enter the 1st sequence x(n)=')
h=input('enter the 1st sequence h(n)=')
N=length(x)+length(h)-1;

%zero padding the sequences x & h

x=[x zeros(1,N-length(x))]
h=[h zeros(1,N-length(h))]
%computing linear convolution using circular convolution equation
y=cconv(x,h,N); %circular convolution of x & h
disp('circular convolution of x & h is y=');
disp(y);

disp(‘Verification of the result’)

Dept. of EEE, Citech Page 22

 Digital signal processing laboratory (15EEL68)

z=conv(x,h)

n=0:N-1;
stem(n,y);
title('Circular convolution output y(n)');

Result:
Enter the sequence 1: [1 2 ]
Enter the sequence 2: [1 3 4 6 7]

circular convolution of x & h is y=

1 5 10 14 19 14
Verification of the result
1 5 10 14 19 14

Solved examples:
Example.1:
Let’s take x1(n) = {1, 1, 2, 1}, x2(n) = {1, 2, 3, 4} and N
Arrange x1(n) and x2(n) in circular fashion as shown below:

To get x2(-m), rotate x2(m) by 4 samples in clockwise direction.

Dept. of EEE, Citech Page 23

 Digital signal processing laboratory (15EEL68)

x3(0) = x1(m) x2((1-m))N

= x1(0) x2(0) + x1(1) x2(3) + x1(2) x2(2) + x1(3) x2(1) = 1 + 4 + 6 + 2

x3(0) = 13

Keep x1(m) constant and rotate x2(-m) once to compute further values.
To get x3(1) rotate x2(-m) by one sample in anti-clockwise direction

x2(1-m)

x3(1) = x1(m) x2((1-m))N

= x1(0) x2(1) + x1(1) x2(0) + x1(2) x2(3) + x1(3) x2(2)= 2 + 1 + 8 + 3
x3(1) = 14

To get x3(2) rotate x2(1-m) by one sample in anti-clockwise direction x2(2-m)

x3(2) = x1(m) x2((2-m))N

= x1(0) x2(2) + x1(1) x2(1) + x1(2) x2(0) + x1(3) x2(3) = 3 + 2 + 2 + 4
x3(1) = 11
To get x3(3) rotate x2(2-m) by one sample in anti-clockwise direction
x2(3-m)

Dept. of EEE, Citech Page 24

 Digital signal processing laboratory (15EEL68)

x3(3) = x1(m) x2((2-m))N

= x1(0) x2(3) + x1(1) x2(2) + x1(2) x2(1) + x1(3) x2(0)=4 + 3 + 4+1
x3(3) = 12

The convoluted signal is,

x3(n) = {13, 14, 11, 12}

Output:
Enter the first sequence: [1 1 2 1]
Enter the second sequence: [1 2 3 4]
The circular convolution of two given sequences is:
13 14 11 12
Matrix method: x(n)={1,2,3,4}h(n)={1 ,2 ,3,4 }

y (n)={26,28,26,20}

Viva Questions:

1. What is the mathematical model for circular convolution?

2. Distinguish between circular convolution and linear convolution of two sequences?
3. Give the steps to get the result of linear convolution using the method of circular
convolution?
4. What is zero padding?
5. Define circular convolution
6. What is length of circular convoluted sequence?
7. Explain the following keywords max and cconv.
8. Different methods of calculating circular convolution
9. Difference between linear convolution and circular convolution
10. Write the mathematical equation for the circular convolution

Dept. of EEE, Citech Page 25

 Digital signal processing laboratory (15EEL68)

EXPERIMENT NO 5. Computation of N point DFT of a given sequence and to plot

magnitude and phase spectrum.

Aim: To compute DFT of the given sequence.

Theory:
• Discrete Fourier Transform (DFT) is used for performing frequency analysis of
discrete time signals. DFT gives a discrete frequency domain representation whereas
the other transforms are continuous in frequency domain.
• The N point DFT of discrete time signal x[n] is given by the equation
N -1 − j 2kn

X (k ) =  x[n]e N
; k = 0,1,2,....N - 1
n =0
Where N is chosen such that N  L , where L=length of x[n].
• The inverse DFT allows us to recover the sequence x[n] from the frequency samples.
j 2kn
1 N −1
x[n] =  x[n]e N ; n = 0,1,2,....N - 1
N k =0
• X(k) is a complex number (remember ejw=cosw + jsinw). It has both magnitude and
phase which are plotted versus k. These plots are magnitude and phase spectrum of
x[n]. The ‘k’ gives us the frequency information.
• Here k=N in the frequency domain corresponds to sampling frequency (fs). Increasing
N, increases the frequency resolution, i.e., it improves the spectral characteristics of the
sequence. For example if fs=8kHz and N=8 point DFT, then in the resulting spectrum,
k=1 corresponds to 1kHz frequency. For the same fs and x[n], if N=80 point DFT is
computed, then in the resulting spectrum, k=1 corresponds to 100Hz frequency. Hence,
the resolution in frequency is increased.
• Since N  L , increasing N to 8 from 80 for the same x[n] implies x[n] is still the same
sequence (<8), the rest of x[n] is padded with zeros. This implies that there is no further
information in time domain, but the resulting spectrum has higher frequency resolution.
This spectrum is known as ‘high density spectrum’ (resulting from zero padding x[n]).
Instead of zero padding, for higher N, if more number of points of x[n] are taken (more
data in time domain), then the resulting spectrum is called a ‘high resolution spectrum’.

Algorithm:
1. Input the sequence for which DFT is to be computed.
2. Input the length of the DFT required (say 4, 8, >length of the sequence).
3. Compute the DFT using the ‘fft’ command.
4. Plot the magnitude & phase spectra.

MATLAB Implementation:
The Discrete Fourier transform is calculated by formula by using for loop.The magnitude
spectrum is computed using the function ABS,ABS(X) is the absolute value of the elements of
X. When X is complex, ABS(X) is the complex modulus (magnitude) of the elements of X.The
phase spectrum is computed using the function ANGLE Phase angle. ANGLE (H) returns the
phase angles, in radians, of a matrix with complex elements.
MATLAB:
clc
x=input('enter the time sequence');

Dept. of EEE, Citech Page 26

 Digital signal processing laboratory (15EEL68)

N=input('enter the no.of points of DFT');

for k=0:N-1
X(k+1)=0;
for n=0:length(x)-1
X(k+1)= X(k+1) + x(n+1)* exp (-j*2*pi*n*k/N)
end
end
n=0:length(x)-1;
stem(n,x);
title('time sequence');
xlabel('n');
ylabel('x(n)');

k=0:N-1;
subplot(2,1,1);
stem(k,abs(X));
title('magnitude response');
xlabel('k');
ylabel('mag(X(k))');

subplot(2,1,2);
stem(k,angle(X));
title('phase response');
xlabel('k');
ylabel('phase(X(k))');

RESULT:

enter the time sequence [1 2 3 4]

enter the no. of points of DFT 4

10.0000 + 0.0000i
-2.0000 + 2.0000i
-2.0000 - 0.0000i
-2.0000 - 2.0000i

Dept. of EEE, Citech Page 27

 Digital signal processing laboratory (15EEL68)

Solved examples:

Example:
Let us assume the input sequence x[n] = [1 1 0 0]

We have,

k = 0, 1…………, N-1

For k = 0, N=4 = x(0) + x(1) + x(2) + x(3)

X(0) = 1+1+0+0 = 2

For k = 1 =x(0) + x(1) e-jπ/2 + x(2) e-jπ + x(3) e-j3π/2

= 1 + cos(π/2) - jsin(π/2)
X(1) = 1 – j

For k = 2, = x(0) + x(1) e-jπ + x(2) e-j2π + x(3) e-j3π

= 1 + cos π – jsin π
X(2) = 1-1 = 0

For k = 3, = x(0) + x(1) e-j3π/2 + x(2) e-j3π + x(3) e-j9π/2

= 1 + cos(3π/2) - jsin(3π/2)
X(3) = 1 + j

The DFT of the given sequence is,

X(k) = { 2, 1-j, 0, 1+j }

Dept. of EEE, Citech Page 28

 Digital signal processing laboratory (15EEL68)

To find Magnitude of X(k):

Magnitude= (a2+b2)1/2
Where a and b are real and imaginary parts respectively
Magnitude of X(k)=2.0000 1.4142 0 1.4142

To fine Phase of X(k):

Phase=tan-1(b/a)
Phase of X(k)=0 -0.7854 0 0.7854

Viva questions:
1. Mathematical representations of Fourier transform of sequence.
2. Mathematical representations of discrete Fourier transform of sequence.
3. What is spectrum ?
4. Why zero padding necessary ?
5. What is frequency resolution ?
6. How to calculate magnitude and phase values of a complex number ?
7. Explain following keywords abs, angle and cos.

Dept. of EEE, Citech Page 29

 Digital signal processing laboratory (15EEL68)

EXPERIMENT NO 6: Linear and circular convolution of two given sequences using

DFT and IDFT.

Aim: To perform linear convolution and circular convolution of two given sequences x1 & x2
using DFT & IDFT

Theory:
• By multiplying two N-point DFTs in the frequency domain, we get the circular
convolution in the time domain.
N − DFT
y[n] = x[n]  h[n] ⎯ ⎯⎯→ Y (k ) = X (k )  H (k )
• Circular Convolution is made identical to linear convolution by choosing the number
of points in DFT as N ≥ xlength + hlength -1.
• For Circular Convolution, N = max (xlength, hlength).

Algorithm:
1. Input the two sequences as xn1, xn2
2. Compute N = xn1 length + xn2 length -1
3. Obtain N-point DFTs (Xf1, Xf2) of both the sequences.
4. Multiply both the DFTs (Y=Xf1*Xf2).
5. Perform IDFT on Y to get y[n] (the linearly convolved sequence) in time domain.
6. Verify using the ‘conv’ command.
7. Plot the sequences

MATLAB Implementation:
MATLAB has the function fft(x,N) to perform the N point DFT. The function IFFT(X)
performs the inverse discrete Fourier transform of X. IFFT(X,N) is the N-point inverse
transform.
MAX - Largest component. For vectors, MAX(X) is the largest element in X. LENGTH
Length of a vector. LENGTH(X) returns the length of vector X. In MATLAB, the operator *
implies Matrix multiplication. X*Y is the matrix product of X and Y. The number of columns
of X must be equal to the number of rows of Y.The operator .* Array multiplication. X.*Y
denotes element-by-element multiplication. X and Y must have the same dimensions (unless
one is a scalar). In the below program to multiply the two DFTs (vectors), we use .* (dot star)
operator.

Dept. of EEE, Citech Page 30

 Digital signal processing laboratory (15EEL68)

MATLAB Program:
%Linear Convolution Program using DFT
x = input('enter the first sequence');
h= input('enter the second sequence');
N=length(x)+length(h)-1;
xk= fft(x,N); %obtain DFTs
hk= fft(h,N);
yk= xk.*hk; %Perform vector multiplication
yn= ifft(yk,N); %ifft to get y[n]
disp('linear convolution of x & h is yn=');
disp(yn);
stem(yn);
xlabel(' n');
ylabel('y[n]');
title('linear convolution output y(n)');

Result:
enter the first sequence [1 2 3 4]
enter the second sequence [1 2 3 4]

linear convolution of x & h is yn=

1.0000 4.0000 10.0000 20.0000 25.0000 24.0000 16.0000
output using conv=
1 4 10 20 25 24 16
The output plot is shown in Figure10.1 window
linear convolution output y(n)
25

20

15
y[n]

10

0
1 2 3 4 5 6 7
n

Figure10.1 Linear convolution output y(n)

Dept. of EEE, Citech Page 31

 Digital signal processing laboratory (15EEL68)

MATLAB Program:
%Circular Convolution Program using DFT
x = input('enter the first sequence');
h= input('enter the second sequence');
N=max(length(x),length(h));
xk= fft(x,N); %obtain DFTs
hk= fft(h,N);
yk= xk.*hk; %Perform vector multiplication
yn= ifft(yk,N); %ifft to get y[n]
disp('circular convolution of x & h is yn=');
disp(yn);
stem(yn);
xlabel(' n');
ylabel('y[n]');
title('linear convolution output y(n)');

Result:
enter the first sequence [1 2 3 4]
enter the second sequence [1 2 3 4]

circular convolution of x & his yn=

26 28 26 20

The output plot is shown in Figure11.1 window

Circular convolution output y(n)

30

25

20
y[n]

15

10

0
1 1.5 2 2.5 3 3.5 4
n

Figure11.1 Circular convolution output

Solved examples for linear convolution

Example.1:
x1(n) = {1, 1, 2}
x2(n) = {1, 2}

For linear convolution,

Length N = Length(x1) + Length(x2) - 1

Dept. of EEE, Citech Page 32

 Digital signal processing laboratory (15EEL68)

N=3+2–1=4

Convolution of two sequences x1(n) and x2(n) is,

x3(n) = IDFT[X3(k)]
x3(n) = IDFT[X1(k) X2(k)]

Where,
X1(k) = DFT [x1(n)]
X2(k) = DFT [x2(n)]

k= 0, 1, 2, … , N-1

Given x1(n) = {1, 1, 2} and N=4

=1+1+2=4

= 1 – j – 2 = -1 – j

n
=1–1+2=2

X1(3) = 1 + j – 2 = -1 + j

X1(k) = {4, -1-j, 2, -1+j}

Now,
k= 0, 1, 2, … , N-1

Given x2(n) = {1, 2} and N=4

=1+2=3

= 1 + 2(-j) = 1 - j2

= 1 + 2(-1) = -1

/2
= 1 + 2(j) = 1 + j2

X2(k) = {3, 1-j2, -1, 1+j2}

We know that,

X3(k) = X1(k) X2(k)

X3(k) = {12, -3+j, -2, -3-j}

Convolution of two given sequences is,

Dept. of EEE, Citech Page 33

 Digital signal processing laboratory (15EEL68)

x3(n) = IDFT[X3(k)]

n = 0, 1, 2, …. , N-1

= (1/4) [12 – 3 +j -2 -3 –j] = 1

= (1/4) [12 + (-3 + j) j + (-2) (-1) + (-3 - j) (-j)] = 3

= (1/4) [12 + (-3 + j) (-1) + (-2) (1) + (-3 - j) (-1)] = 4

4=
(1/4) [12 + (-3 + j) (-j) + (-2) (-1) + (-3 - j) (j)] = 4

Convoluted sequence of two given sequences is,

x3(n) = {1, 3, 4, 4}

Output:
Enter the 1st seq: [1 1 2]
Enter the 2nd seq: [1 2]
Linear convolution using DFT and IDFT:
1 3 4 4

Solved examples for circular convolution :

Circular convolution:
Let x1(n) and x2(n) are finite duration sequences both of length N with DFT’s X1(k) and
X2(k). Convolution of two given sequences x1(n) and x2(n) is given by,
x3(n) = IDFT[X3(k)]
x3(n) = IDFT[X1(k) X2(k)]

Where,
X1(k) = DFT [x1(n)]
X2(k) = DFT [x2(n)]

Example:
x1(n) = {1, 1, 2, 1}
x2(n) = {1, 2, 3, 4}

X1(k) k= 0, 1, 2, … , N-1

Given x1(n) = {1, 1, 2, 1} and N=4

=1+1+2+1=5

= 1 – j – 2 + j = -1

=1–1+2-1=1

Dept. of EEE, Citech Page 34

 Digital signal processing laboratory (15EEL68)

= 1 + j - 2 - j = -1

X1(k) = {5, -1, 1, -1}

Now, k= 0, 1, 2, … , N-1

Given x2(n) = {1, 2, 3, 4} and N=4

= 1 + 2 + 3 + 4 = 10

= 1 + 2(-j) + 3(-1) + 4(j) = -2 + j2

= 1 + 2(-1) + 3(1) + 4(-1) = -2

= 1 + 2(j) + 3(-1) + 4(-j) = -2 – j2

X2(k) = {10, -2+j2, -2, -2-j2}

We know that,

X3(k) = X1(k) X2(k)

X3(k) = {50, 2 – j2, -2, 2 + j2}

Convolution of two given sequences is,

x3(n) = IDFT[X3(k)]

n = 0, 1, 2, …. , N-1

= (1/4) [50 + 2 - j2 – 2 + 2 + j2] = 13

= (1/4) [50 + (2 - j2) j + (-2) (-1) + (2 + j2) (-j)] = 14

= (1/4) [50 + (2 - j2) (-1) + (-2) (1) + (2 + j2) (-1)] = 11

4=
(1/4) [50 + (2 - j2) (-j) + (-2) (-1) + (2 + j2) (j)] = 12

Convoluted sequence of two given sequences is,

x3(n) = {13, 14, 11, 12}

Dept. of EEE, Citech Page 35

 Digital signal processing laboratory (15EEL68)

Viva questions:

1. Explain following keywords fft and ifft

2. Different methods of calculating circular convolution
3. Define circular convolution
4. Functions use to calculate circular convolution
5. What is length of circular convoluted sequence?

Dept. of EEE, Citech Page 36

 Digital signal processing laboratory (15EEL68)

EXPERIMENT NO 7: Solution of a difference equation.

Aim: To obtain the impulse response/step response of a system described by the given
difference equation.
Theory:
• A difference equation with constant coefficients describes a LTI system. For example
the difference equation y[n] + 0.8y[n-2] + 0.6y[n-3] = x[n] + 0.7x[n-1] + 0.5x[n-2]
describes a LTI system of order 3. The coefficients 0.8, 0.7, etc are all constant i.e.,
they are not functions of time (n). The difference equation y[n]+0.3ny[n-1]=x[n]
describes a time varying system as the coefficient 0.3n is not constant.
• The difference equation can be solved to obtain y[n], the output for a given input x[n]
by rearranging as y[n] = x[n] + 0.7x[n-1]+0.5x[n-2]- 0.8y[n-2]- 0.6y[n-3] and solving.
• The output depends on the input x[n]
o With x[n]= δ[n], an impulse, the computed output y[n] is the impulse
response.
o If x[n]=u[n], a step response is obtained.
o If x[n] = cos(wn) is a sinusoidal sequence, a steady state response is obtained
(wherein y[n] is of the same frequency as x[n], with only an amplitude gain and
phase shift-refer Fig.7.3).
o Similarly for any arbitrary sequence of x[n], the corresponding output response
y[n] is computed.
• The difference equation containing past samples of output, i.e., y[n-1], y[n-2], etc leads
to a recursive system, whose impulse response is of infinite duration (IIR). For such
systems the impulse response is computed for a large value of n, say n=100 (to
approximate n=∞). The MATLAB function filter is used to compute the impulse
response/ step response/ response to any given x[n]. Note: The filter function evaluates
the convolution of an infinite sequence (IIR) and x[n], which is not possible with conv
function (remember conv requires both the sequences to be finite).
• The difference equation having only y[n] and present and past samples of input (x[n],
x[n-k]), represents a system whose impulse response is of finite duration (FIR). The
response of FIR systems can be obtained by both the ‘conv’ and ‘filter’ functions. The
filter function results in a response whose length is equal to that of the input x[n],
whereas the output sequence from conv function is of a longer length (xlength +
hlength-1).

Algorithm:
1. Input the two sequences as a and b representing the coefficients of y and x.
2. If IIR response, then input the length of the response required (say 100, which can be
made constant).
3. Compute the output response using the ‘filter’ command.
4. Plot the input sequence & impulse response (and also step response, etc if required).

Dept. of EEE, Citech Page 37

 Digital signal processing laboratory (15EEL68)

MATLAB Implementation:
MATLAB has an inbuilt function ‘filter’ to solve difference equations numerically, given
the input and difference equation coefficients (b,a).
y=filter(b,a,x)
where x is the input sequence, y is the output sequence which is of same length as x.
Given a difference equation a0y[n]+a1y[n-1]+a2y[n-2]=b0x[n]+b2x[n-2], the coefficients are
written in a vector as b=[b0 0 b2] and a=[a0 a1 a2]. Note the zero in b (x[n-1] term is missing).
Also remember a0, the coefficient of y[n] should always be 1.
For impulse responsethe numb x[n] = {1,0,0,0,....} er of zeros = the length of the IIR response


required (say 100 implemented by function zeros(1,99)).

For step response x[n] = {1,1,1,1,1,....} the number of ones = the length of the IIR response


required-1 (say 100 implemented by function ones(1,100).

Similarly for any given x[n] sequence, the response y[n] can be calculated.
Given Problem
1) Difference equation y(n) - y(n-1) + 0.9y(n-2) = x(n);

Calculate impulse response h(n) and also step response at n=0,…..,100

2 ) Plot the steady state response y[n] to x[n]=cos(0.05πn)u(n), given y[n]-0.8y[n-1]=x[n]

%To find Impulse Response

N=input('Length of response required=');
b=[1]; %x[n] coefficient
a=[1,-1,0.9]; %y coefficients
x=[1,zeros(1,N-1)]; %impulse input
n=0:1:N-1; %time vector for plotting
h=filter(b,a,x); %impulse response
subplot(2,1,1);
stem(n,x);
title('impulse input');
xlabel('n');
ylabel('δ(n)');
subplot(2,1,2);
stem(n,h);
title('impulse response');
xlabel('n');
ylabel('h(n)');

Result:
Length=100
Plots in Fig. 7.1

Dept. of EEE, Citech Page 38

 Digital signal processing laboratory (15EEL68)

Fig.7.1 Impulse Response

%To find step response

N=input('Length of response required=');
b=[1]; %x[n] coefficient
a=[1,-1,0.9]; %y coefficients
x=[ones(1,N)]; %step input
n=0:1:N-1; %time vector for plotting
y=filter(b,a,x); %step response
subplot(2,1,1);
stem(n,x);
title('step input');
xlabel('n');
ylabel('u(n)');
subplot(2,1,2);
stem(n,h);
title('step response');
xlabel('n');
ylabel('y(n)');

Dept. of EEE, Citech Page 39

 Digital signal processing laboratory (15EEL68)

Result:
Length=100
Plots in Fig. 7.2

Fig.7.2 Step Response

%To find steady state response

%y[n]-0.8y[n-1]=x[n]
N=input('Length of response required=');
b=[1]; %x[n] coefficient
a=[1,-0.8]; %y coefficients
n=0:1:N-1; %time vector
x=cos(0.05*pi*n); % sinusoidal input
y=filter(b,a,x); %steady state response
subplot(2,1,1);
stem(n,x);
title('sinusodial input');
xlabel('n');
ylabel('x(n)');
subplot(2,1,2);
stem(n,h);
title('steady state response');
xlabel('n');
ylabel('y(n)');

Dept. of EEE, Citech Page 40

 Digital signal processing laboratory (15EEL68)

Result:
Length=100
Plots in Fig. 7.3

Fig. 7.3 Steady state Response to a cosine input

Solved examples:

A General Nth order Difference equations looks like, y[n] + a1y[n-1] + a2y[n-2] + …+
aNy[n-N) = b0x[n] + b1x[n-1] + …. + bMx[n-M] Here y[n] is “Most advanced” output
sample and y[n-m] is “Least advanced” output sample. The difference between these two
index values is the order of the difference equation.
Here we have: n-(n-N) = N
We can rewrite the above equation as,
y[n] + Σa[i] y[n-i] = Σbi x[n-i] y[n] becomes,
y[n] = -Σa[i] y[n-i] + Σbi x[n-i]
Example:
y[n+2] – 1.5y[n+1] + y[n] = 2x[n]
In general we start with the “Most advanced” output sample. Here it is y[n+2].
So, here we need to subtract 2 from each sample argument. We get
y[n] – 1.5y[n-1] + y[n-2] = x[n]

Dept. of EEE, Citech Page 41

 Digital signal processing laboratory (15EEL68)

y[n] = 1.5y[n-1] - y[n-2] + x[n]

Let us assume our input x[n] = u[n] = 0 x<0
1 x≥0

In our example we have taken x[n] = u[n] = 0 x<0

1 0≤x<10
We need N past values to solve Nth order difference equation.
y[-2] = 2
y[-1] = 1
Compute y[n] for various values of n
y[0] = 1.5y[-1] - y[-2] + x[0]
= 1.5*1 -2+1
y[0] = 0.5
y[1] = 1.5y[0] - y[-1] + x[1]
= 1.5*0.5 - 1 + 1
y[1] =0.75
y[2] = 1.5y[1] - y[0] + x[2]
= 1.5*0.75 – 0.5 + 1
y[3] = 1.625 y[4] = 2.6875 y[5] = 3.4063 y[6] = 3.4219
similarly compute for n=7,8,9,etc……

Output:
Solution for the given difference equation:
Columns 1 through 7
0.5000 0.7500 1.6250 2.6875 3.4063 3.4219 2.7266
Columns 8 through 10
1.6680 0.7754 0.4951

Viva questions:
1. What is the use of difference equation
2. Define step response
3. Define impulse response
4. Define steady state response
5. Explain following keywords zeros, ones and filter
6. What is the difference between difference equation and differential equation?
7. Differentiate between conv and filter
8. What is meant by zero response?
9. What is meant by forced response?
10. What is transfer function of system?
11. Write the difference equation for an FIR system
12. Write the difference equation for IIR system

Dept. of EEE, Citech Page 42

 Digital signal processing laboratory (15EEL68)

MATLAB Program:
Clc;
x= input('enter the time sequence');
N= input ('enter the no. of points of DFT');
X=fft(x,N)

x1 =ifft(X,N)
n=0:length(x)-1;
stem(n,x);
title('time sequence');
xlabel('n');

ylabel('x(n)');
k=0:N-1;
figure
subplot(2,1,1);
stem(k,abs(X));

title('magnitude response');
xlabel('k');
ylabel('mag(X(k))');

subplot(2,1,2);

stem(k,angle(X));
title('phase response');
xlabel('k');
ylabel('phase(X(k))');

RESULT
enter the time sequence [1 2 2 1 2 1 1 2]
enter the no. of points of DFT8
X=
12 0.4142 – j 0 -2.4142 + j

Dept. of EEE, Citech Page 43

 Digital signal processing laboratory (15EEL68)

0 -2.4142 - j 0 0.4142 + j
x1 =

1 2 2 1 2 1 1 2

Dept. of EEE, Citech Page 44

 Digital signal processing laboratory (15EEL68)

EXPERIMENT NO 9: Design and implementation of IIR filter to meet given

specifications.

Aim: To design and implement an IIR filter for given specifications.

Theory:
There are two methods of stating the specifications as illustrated in previous program. In the
first program, the given specifications are directly converted to digital form and the designed
filter is also implemented. In the last two programs the butterworth and chebyshev filters are
designed using bilinear transformation (for theory verification).
Method I: Given the order N, cutoff frequency fc, sampling frequency fs and the IIR filter type
(butterworth, cheby1, cheby2).
• Step 1: Compute the digital cut-off frequency Wc (in the range -π < Wc < π, with π
corresponding to fs/2) for fc and fs in Hz. For example let fc=400Hz, fs=8000Hz
Wc = 2*π* fc / fs = 2* π * 400/8000 = 0.1* π radians
For MATLAB the Normalized cut-off frequency is in the range 0 and 1, where 1
corresponds to fs/2 (i.e.,fmax)). Hence to use the MATLAB commands
wc = fc / (fs/2) = 400/(8000/2) = 0.1
Note: if the cut off frequency is in radians then the normalized frequency is computed as
wc = Wc / π
• Step 2: Compute the Impulse Response [b,a] coefficients of the required IIR filter and the
response type (lowpass, bandpass, etc) using the appropriate butter, cheby1, cheby2
command. For example given a butterworth filter, order N=2, and a high pass response, the
coefficients [b,a] of the filter are computed using the MATLAB inbuilt command ‘butter’
as [b,a] =butter(N, wc , 'high');

Method 2:
Given the pass band (Wp in radians) and Stop band edge (Ws in radians) frequencies, Pass
band ripple Rp and stopband attenuation As.
• Step 1: Since the frequencies are in radians divide by π to obtain normalized frequencies
to get wp=Wp/pi and ws=Ws/pi
If the frequencies are in Hz (note: in this case the sampling frequency should be given),
then obtain normalized frequencies as wp=fp/(fs/2), ws=fstop/(fs/2), where fp, fstop
and fs are the passband, stop band and sampling frequencies in Hz
• Step 2: Compute the order and cut off frequency as
[N, wc] = BUTTORD(wp, ws, Rp, Rs)
• Step 3: Compute the Impulse Response [b,a] coefficients of the required IIR filter and
the response type as [b,a] =butter(N, wc , 'high');

IMPLEMENTATION OF THE IIR FILTER

1. Once the coefficients of the IIR filter [b,a] are obtained, the next step is to simulate an
input sequence x[n], say input of 100, 200 & 400 Hz (with sampling frequency of fs),
each of 20/30 points. Choose the frequencies such that they are >, < and = to fc.
2. Filter the input sequence x[n] with Impulse Response, to obtain the output of the filter
y[n] using the ‘filter’ command.
3. Infer the working of the filter (low pass/ high pass, etc).

Dept. of EEE, Citech Page 45

 Digital signal processing laboratory (15EEL68)

MATLAB IMPLEMENTATION

BUTTORD Butterworth filter order selection.

[N, Wn] = BUTTORD(Wp, Ws, Rp, Rs) returns the order N of the lowest order digital
Butterworth filter that loses no more than Rp dB in the passband and has at least Rs dB of
attenuation in the stopband. Wp and Ws are the passband and stopband edge frequencies,
normalized from 0 to 1 (where 1 corresponds to pi radians/sample). For example,
Lowpass: Wp = .1, Ws = .2 Highpass: Wp = .2, Ws = .1
Bandpass: Wp = [.2 .7], Ws = [.1 .8] Bandstop: Wp = [.1 .8], Ws = [.2 .7]
BUTTORD also returns Wn, the Butterworth natural frequency (or, the "3 dB frequency")
to use with BUTTER to achieve the specifications.
[N, Wn] = BUTTORD(Wp, Ws, Rp, Rs, 's') does the computation for an analog filter, in
which case Wp and Ws are in radians/second.
When Rp is chosen as 3 dB, the Wn in BUTTER is equal to Wp in BUTTORD.
BUTTER Butterworth digital and analog filter design.
[B,A] = BUTTER(N,Wn) designs an Nth order lowpass digital Butterworth filter and
returns the filter coefficients in length N+1 vectors B (numerator) and A (denominator). The
coefficients are listed in descending powers of z. The cutoff frequency Wn must be 0.0 < Wn
< 1.0, with 1.0 corresponding to half the sample rate. If Wn is a two-element vector, Wn =
[W1 W2], BUTTER returns an order 2N bandpass filter with passband W1 < W < W2.
[B,A] = BUTTER(N,Wn,'high') designs a highpass filter. [B,A] = BUTTER(N,Wn,'stop') is a
bandstop filter if Wn = [W1 W2].
BUTTER(N,Wn,'s'), BUTTER(N,Wn,'high','s') and BUTTER(N,Wn,'stop','s') design analog
Butterworth filters. In this case, Wn is in [rad/s] and it can be greater than 1.0.

Program (Design & implementation)

%generate filter coefficients for the given %order & cutoff Say N=2, fc=150Hz,
%fs=1000 Hz, butterworth filter

[b,a]=butter(2, 150/(1000/2));

%generate simulated input of 100, 300 & 170 Hz, each of 30 points
n=1:30;
f1=100;f2=300;f3=170;fs=1000;
x=[];
x1=sin(2*pi*n*f1/fs);
x2=sin(2*pi*n*f2/fs);
x3=sin(2*pi*n*f3/fs);
x=[x1 x2 x3];
subplot(2,1,1);
stem(x);
title('input');
%generate o/p
y=filter(b,a,x);
subplot(2,1,2);
stem(y);
title('output');
Result:

Dept. of EEE, Citech Page 46

 Digital signal processing laboratory (15EEL68)

Plot is in Fig. 9.1, which shows that 100 Hz is passed, while 300 is cutoff and 170 has slight
attenuation.

Note: If fp,fstp,fs, rp,As are given then use

[N,wc]=buttord(2*fp/fs,2*fstp/fs,rp,As)
[b,a]=butter(N,wc);
If wp & ws are in radians
[N,wc]=buttord(wp/pi,ws/pi,rp,As)
[b,a]=butter(N,wc);
If wc is in radians & N is given
[b,a]=butter(N,wc/pi);
For a bandpass output
wc=[150/(1000/2) 250/(1000/2)];
[b,a]=butter(4, wc);
Plot is in Fig.9.2 where only 170 is passed, while 100 & 300 are cutoff.

Fig 9.1. Low pass IIR filter (butterworth) output

Dept. of EEE, Citech Page 47

 Digital signal processing laboratory (15EEL68)

Fig.9.2 IIR output for bandpass (150-250 Hz)

Viva Questions:

1) What are IIR filters? What does "IIR" mean?

2) Why is the impulse response "infinite"?

3) What is the alternative to IIR filters?

4) What are the advantages of IIR filters (compared to FIR filters)?

5) What are the disadvantages of IIR filters (compared to FIR filters)?

6) What is prewarping?

7) Give the frequency transformation techniques?

8) Give the disadvantage of impulse invariance method?

9) What is an elliptic filter?

10) Differentiate the three protoype filters

a) Based on frequency characteristics b) Based on poles and zeros

11) What is the condition for stability of an IIR system?

12) So where did the infinite impulse response go in the IIR filter equation...........

Dept. of EEE, Citech Page 48

 Digital signal processing laboratory (15EEL68)

Programs for designing of IIR filters (for theory practice):

% Butterworth filter
Given data: rp=1, rs=40, w1=800, w2=1200,ws=3600;
% Analog frequency
aw1= 2*pi*w1/ws;
aw2=2*pi*w2/ws;
% Prewrapped frequency
pw1 = 2*tan(aw1/2);
pw2 = 2*tan(aw2/2);
%Calculate order and cutoff freq
[n,wc]= buttord (pw1,pw2,rp,rs,'s');
% analog filter transfer
[b,a] = butter(n,wc,'s');
% obtaining the digital filter using bilinear transformation
fs=1;
[num,den]= bilinear(b,a,fs);
%plot the frequency response
[mag,freq1]=freqz(num,den,128);
freq=freq1*ws/(2*pi);
m = 20*log10(abs(mag));
plot(freq,m);
grid;
Result:
rp = 1
rs = 40
w1 = 800
w2 = 1200
plot is in Fig. 9.3
To design a chebyshev filter for given specifications
%Given data
rp=1,rs=40,w1=800,w2=1200,ws=3600

%Analog frequencies
aw1= 2*pi*w1/ws;
aw2=2*pi*w2/ws;

% Prewrapped frequency assuming T=1/fs

pw1 = 2*tan(aw1/2);
pw2 = 2*tan(aw2/2);

[n,wc]= cheb1ord (pw1,pw2,rp,rs,'s');

[b,a] = cheby1(n,rp,wc,'s');

% obtaining the digital filter using bilinear transformation

fs=1;
[num,den]= bilinear(b,a,fs);

%plot the frequency response

[mag,freq1]=freqz(num,den,128);

Dept. of EEE, Citech Page 49

 Digital signal processing laboratory (15EEL68)

freq=freq1*ws/(2*pi);
m = 20*log10(abs(mag));
plot(freq,m);
grid;
Result:
rp = 1, rs = 40
w1 = 800, w2 = 1200
ws = 3600
Plot is in Fig 9.4

Fig 9.3. Butterworth Filter

Fig.9.4 Chebyshev Filter

Dept. of EEE, Citech Page 50

 Digital signal processing laboratory (15EEL68)

EXPERIMENT NO 10: Design and implementation of FIR filter to meet given

specifications.
Aim: To design and implement a FIR filter for given specifications.
Theory:
There are two types of systems – Digital filters (perform signal filtering in time domain) and
spectrum analyzers (provide signal representation in the frequency domain). The design of a
digital filter is carried out in 3 steps- specifications, approximations and implementation.

DESIGNING AN FIR FILTER (using window method):

Method I: Given the order N, cutoff frequency fc, sampling frequency fs and the window.
• Step 1: Compute the digital cut-off frequency Wc (in the range -π < Wc < π, with π
corresponding to fs/2) for fc and fs in Hz. For example let fc=400Hz, fs=8000Hz
Wc = 2*π* fc / fs = 2* π * 400/8000 = 0.1* π radians
For MATLAB the Normalized cut-off frequency is in the range 0 and 1, where 1
corresponds to fs/2 (i.e.,fmax)). Hence to use the MATLAB commands
wc = fc / (fs/2) = 400/(8000/2) = 0.1
Note: if the cut off frequency is in radians then the normalized frequency is computed
as wc = Wc / π
• Step 2: Compute the Impulse Response h(n) of the required FIR filter using the given
Window type and the response type (lowpass, bandpass, etc). For example given a
rectangular window, order N=20, and a high pass response, the coefficients (i.e., h[n]
samples) of the filter are computed using the MATLAB inbuilt command ‘fir1’ as
h =fir1(N, wc , 'high', boxcar(N+1));
Note: In theory we would have calculated h[n]=hd[n]×w[n], where hd[n] is the desired
impulse response (low pass/ high pass,etc given by the sinc function) and w[n] is the
window coefficients. We can also plot the window shape as stem(boxcar(N)).
Plot the frequency response of the designed filter h(n) using the freqz function and observe the
type of response (lowpass / highpass /bandpass).

Method 2:
Given the pass band (wp in radians) and Stop band edge (ws in radians) frequencies, Pass
band ripple Rp and stopband attenuation As.
• Step 1: Select the window depending on the stopband attenuation required. Generally
if As>40 dB, choose Hamming window. (Refer table )
• Step 2: Compute the order N based on the edge frequencies as
Transition bandwidth = tb=ws-wp;
N=ceil (6.6*pi/tb);
• Step 3: Compute the digital cut-off frequency Wc as
Wc=(wp+ws)/2
Now compute the normalized frequency in the range 0 to 1 for MATLAB as
wc=Wc/pi;
Note: In step 2 if frequencies are in Hz, then obtain radian frequencies (for computation of tb
and N) as wp=2*pi*fp/fs, ws=2*pi*fstop/fs, where fp, fstop and fs are the passband, stop band
and sampling frequencies in Hz
• Step 4: Compute the Impulse Response h(n) of the required FIR filter using N, selected
window, type of response(low/high,etc) using ‘fir1’ as in step 2 of method 1.

Dept. of EEE, Citech Page 51

 Digital signal processing laboratory (15EEL68)

IMPLEMENTATION OF THE FIR FILTER

1. Once the coefficients of the FIR filter h[n] are obtained, the next step is to simulate an
input sequence x[n], say input of 100, 200 & 400 Hz (with sampling frequency of fs),
each of 20/30 points. Choose the frequencies such that they are >, < and = to fc.
2. Convolve input sequence x[n] with Impulse Response, i.e., x (n)*h (n) to obtain the
output of the filter y[n]. We can also use the ‘filter’ command.
3. Infer the working of the filter (low pass/ high pass, etc).

MATLAB IMPLEMENTATION

FIR1 Function

B = FIR1(N,Wn) designs an N'th order lowpass FIR digital filter and returns the filter
coefficients in length N+1 vector B. The cut-off frequency Wn must be between 0 < Wn < 1.0,
with 1.0 corresponding to half the sample rate. The filter B is real and has linear phase, i.e.,
even symmetric coefficients obeying B(k) = B(N+2-k), k = 1,2,...,N+1.
If Wn is a two-element vector, Wn = [W1 W2], FIR1 returns an order N bandpass filter with
passband W1 < W < W2. B = FIR1(N,Wn,'high') designs a highpass filter. B =
FIR1(N,Wn,'stop') is a bandstop filter if Wn = [W1 W2]. If Wn is a multi-element vector,
Wn = [W1 W2 W3 W4 W5 ... WN], FIR1 returns an order N multiband filter with bands
0 < W < W1, W1 < W < W2, ..., WN < W < 1.
FREQZ Digital filter frequency response. [H,W] = FREQZ(B,A,N) returns the N-point
complex frequency response vector H and the N-point frequency vector W in radians/sample
of the filter whose numerator and denominator coefficients are in vectors B and A. The
frequency response is evaluated at N points equally spaced around the upper half of the unit
circle. If N isn't specified, it defaults to 512.
For FIR filter enter A=1 and B = h[n] coefficients. Appropriately choose N as 128, 256, etc

Table: Commonly used window function characteristics

___________________________________________________________________________
Window Transition Width  Min. Stop band Matlab
Name Approximate Exact values Attenuation Command
___________________________________________________________________________
4 1.8
Rectangular 21db B = FIR1(N,Wn,boxcar)
M M
8 6.1
Bartlett 25db B = FIR1(N,Wn,bartlett)
M M
8 6.2
Hanning 44db B = FIR1(N,Wn,hanning)
M M
8 6.6
Hamming 53db B= FIR1(N,Wn,hamming)
M M

Dept. of EEE, Citech Page 52

 Digital signal processing laboratory (15EEL68)

12 11
Blackman 74db B = FIR1(N,Wn,blackman)
M M
Program:
%Design and implementation of FIR filter Method 1
%generate filter coefficients for the given %order & cutoff Say N=33, fc=150Hz,
%fs=1000 Hz, Hamming window
h=fir1(33, 150/(1000/2),hamming(34));
%generate simulated input of 50, 300 & 200 Hz, each of 30 points
n=1:30;
f1=50;f2=300;f3=200;fs=1000;
x=[];
x1=sin(2*pi*n*f1/fs);
x2=sin(2*pi*n*f2/fs);
x3=sin(2*pi*n*f3/fs);
x=[x1 x2 x3];
subplot(2,1,1);
stem(x);
title('input');
%generate o/p
%y=conv(h,x);
y=filter(h,1,x);
subplot(2,1,2);
stem(y);
title('output');
Result:
Plots are in Fig.10.1 & 10.2

Notice that frequencies below 150 Hz are passed & above 150 are cutoff
%Method 2: the following program gives only the design of the FIR filter- for
implementation continue with the next program (after h[n])
%input data to be given: Passband & Stopband frequency
% Data given: Passband ripple & stopband attenuation As. If As>40 dB, Choose
hamming
clear
wpa=input('Enter passband edge frequency in Hz');
wsa= input('Enter stopband edge frequency in Hz');
ws1= input('Enter sampling frequency in Hz');
%Calculate transmission BW,Transition band tb,order of the filter
wpd=2*pi*wpa/ws1;
wsd=2*pi*wsa/ws1;
tb=wsd-wpd;
N=ceil(6.6*pi/tb)
wc=(wsd+wpd)/2;
%compute the normalized cut off frequency
wc=wc/pi;
%calculate & plot the window
hw=hamming(N+1);
stem(hw);
title('Fir filter window sequence- hamming window');
%find h(n) using FIR

Dept. of EEE, Citech Page 53

 Digital signal processing laboratory (15EEL68)

h=fir1(N,wc,hamming(N+1));
%plot the frequency response
figure(2);
[m,w]=freqz(h,1,128);
mag=20*log10(abs(m));
plot(ws1*w/(2*pi),mag);
title('Fir filter frequency response');
grid;

Result:
Enter passband edge frequency in Hz100
Enter stopband edge frequency in Hz200
Enter sampling frequency in Hz1000
N = 33
Plots are as in Fig.10.3 and 10.4

Inference:
Notice the maximum stop band attenuation of 53 dB from plot 10.4

Fig.10.1 (using filter command) Fig.10.2 (using conv command)

Fig.10.3 Fig.10.4

Dept. of EEE, Citech Page 54

 Digital signal processing laboratory (15EEL68)

Fig.10.5 Freqs f1=150;f2=300;f3=170;fs=1000;

Dept. of EEE, Citech Page 55

 Digital signal processing laboratory (15EEL68)

EXPERIMENT 11
AIM: To design and test FIR filters based on frequency sampling method
THEORY:
IIR digital filters were designed to give a desired magnitude frequency
response without caring for phase response. In many applications a linear phase is
required throughout the pass band of the filter to preserve the shape of the given
signal in the pass band.
Frequency sampling method consists simply of uniformly sampling the desired
frequency response, and performing an inverse DFT to obtain the corresponding
finite impulse response.
In frequency sampling method we sample the desired frequency response
Hd(w) at N equally sampled points in the interval 0 to 2П where w =( 2Пk/N). This is
clearly a sampling of DTFT to obtain DFT H(k).
H(k) = Hd(w) where k is 0,1…..N-1.
H(k) does not have conjugate symmetry about the folding index N/2. This
symmetry is compulsory for inverse DFT to yield a filter having real coefficients for its
impulse response h(n). Hence we force the conjugate symmetry as follows.
For N= odd, H(0) is real and H(N-k) = H*(k) for k= 1,2,….(N-1)/2
For N=even, H(0)= real and H(N-k)= H*(k) for k= 1,2,….(N/2)-1 and H(N/2)=0
It is suitable in the implementation of parallel processors. It does not require
the computation of inverse DFT which is tedious for long filters. It requires (2N-1)
complex multiplications and one real multiplication per sample whereas FIR needs
only N real multiplications.

RESULT:
LOW PASS

Dept. of EEE, Citech Page 56

 Digital signal processing laboratory (15EEL68)

HIGH PASS

BAND PASS

BAND STOP

Dept. of EEE, Citech Page 57

 Digital signal processing laboratory (15EEL68)

EXPERIMENT 12:

Dept. of EEE, Citech Page 58

 Digital signal processing laboratory (15EEL68)

Cplxpair Sort complex numbers into complex conjugate pairs. Syntax: B = cplxpair(A)
If A is a vector, cplxpair(A) returns A with complex conjugate pairs grouped together.

[r,p,k] = residue(b,a)
[b,a] = residue(r,p,k)
[r,p,k] = residue(b,a) finds the residues, poles, and direct term of a partial fraction expansion
of the ratio

Dept. of EEE, Citech Page 59

 Digital signal processing laboratory (15EEL68)

of two polynomials, and , of the form

Dept. of EEE, Citech Page 60

 Digital signal processing laboratory (15EEL68)

b=[1 0 0 0 16.0625 0 0 0 1]
>> c=roots(b)

c=

-1.4142 + 1.4142i
-1.4142 - 1.4142i
1.4142 + 1.4142i
1.4142 - 1.4142i
-0.3536 + 0.3536i
-0.3536 - 0.3536i
0.3536 + 0.3536i
0.3536 - 0.3536i

>> conv([1 -c(1)], [1 -c(2)])

ans =

1.0000 2.8284 4.0000

>> conv([1 -c(3)], [1 -c(4)])

ans =

1.0000 -2.8284 4.0000

>> conv([1 -c(5)], [1 -c(6)])

ans =

1.0000 0.7071 0.2500

>> conv([1 -c(7)], [1 -c(8)])

ans =

1.0000 -0.7071 0.2500

Dept. of EEE, Citech Page 61