JADAVPUR UNIVERSITY

SALT LAKE CAMPUS, KOLKATA 700 106

Digital Signal Processing Laboratory Manual

Course code: IEE/PC/B/S/321

1
2
Course code: Digital Signal Processing L T P C
IEE/PC/B/S/321 Laboratory 0 0 3 1.5

Course Upon completion of the course, the students will be able to

Outcomes: CO1: Examine and execute MATLAB signal processing functions (S2, A2)
CO2: Examine and execute different mathematical operations on discrete
signals. (S2, A2)
CO3: Examine and execute different digital filters (S2, A2)
CO4: Demonstrate real-time signals and examine their response with
different digital filters (K3,S3-Demonstrate, A2)
List of 1. Find out the impulse and step response of the discrete systems
Experiments 2. Verification of the sampling theorem.
3. Generate different signals using MATLAB
a. Sine wave
b. Cosine wave
c. Unit impulse
d. Unit step wave
e. Square wave
f. Exponential waveform
g. Sawtooth waveform
h. Amplitude-modulated Signal
i. Frequency-modulated signal
4. Study the cross-correlation and auto-correlation of two signals.
5. To find DFT / IDFT and FFTofa given discrete signal
6. Program to obtain Linear and circular Convolution
7. Program for Computing cross-correlation and auto-correlation
8. To find the frequency response of a given system(transfer function/
9. Implementation of FFT of a given sequence
10. Determination of the Power Spectrum of a given signal.
11. Implementation of LP FIR filter for a given sequence
12. Implementation of HP FIR filter for a given sequence
13. Implementation of LP IIR filter for a given sequence
14. Implementation of HP IIR filter for a given sequence
15. Generation of DTMF signals
16. Implementation of the Decimation Process
17. Implementation of Interpolation Process
18. Implementation of I/D sampling rate converters
19. Recovering a sinusoidal signal buried in noise
20. Familiarization with DSP starter kits: Implement an IIR/FIR
filter(LPF/BPF/HPF/ BSF) using a DSK (TMS320C6711).

3
INRODUCTION

MATLAB: MATLAB is a software package for high performance numerical computation and
visualization provides an interactive environment with hundreds of built in functions for
technical computation, graphics and animation. The MATLAB name stands for MATrix
Laboratory

At its core ,MATLAB is essentially a set (a “toolbox”) of routines (called “m files”

or “mex files”) that sit on your computer and a window that allows you to create new variables
with names (e.g. voltage and time) and process those variables with any of those routines (e.g.
plot voltage against time, find the largest voltage, etc).

It also allows you to put a list of your processing requests together in a file and save that
combined list with a name so that you can run all of those commands in the same order at some
later time. Furthermore, it allows you to run such lists of commands such that you pass in data

4
5
and/or get data back out (i.e. the list of commands is like a function in most programming
languages). Once you save a function, it becomes part of your toolbox (i.e. it now looks to you as
if it were part of the basic toolbox that you started with).

For those with computer programming backgrounds: Note that MATLAB runs as an
interpretive language (like the old BASIC). That is, it does not need to be compiled. It simply
reads through each line of the function, executes it, and then goes on to the next line. (In
practice, a form of compilation occurs when you first run a function, so that it can run faster the
next time you run it.)

MATLAB Windows :

MATLAB works with through three basic windows

Command Window : This is the main window .it is characterized by MATLAB command
prompt >> when you launch the application program MATLAB puts you in this window all
commands including those for user-written programs ,are typed in this window at the MATLAB
prompt

Graphics window: the output of all graphics commands typed in the command window are
flushed to the graphics or figure window, a separate gray window with white background color
the user can create as many windows as the system memory will allow

Edit window: This is where you write edit, create and save your own programs in files called M
files.

Input-output:

MATLAB supports interactive computation taking the input from the screen and flushing, the
output to the screen. In addition it can read input files and write output files

Data Type: the fundamental data –type in MATLAB is the array. It encompasses several distinct
data objects- integers, real numbers, matrices, charcter strings, structures and cells.There is no
need to declare variables as real or complex, MATLAB automatically sets the variable to be real.

Dimensioning: Dimensioning is automatic in MATLAB. No dimension statements are required

for vectors or arrays .we can find the dimensions of an existing matrix or a vector with the size
and length commands.

6
Where to work in MATLAB?
All programs and commands can be entered either in the
a) Command window
b) As an M file using Matlab editor
Note: Save all M files in the folder 'work' in the current directory. Otherwise you have to
locate the file during compiling.

Typing quit in the command prompt>> quit, will close MATLAB Matlab Development
Environment.

For any clarification regarding plot etc, which are built in functions type help topic
i.e. help plot

Basic Instructions in Mat lab

1. T = 0: 1:10

This instruction indicates a vector T which as initial value 0 and final value 10 with an
increment of 1

Therefore T = [0 1 2 3 4 5 6 7 8 9 10]

2. F= 20: 1: 100

Therefore F = [20 21 22 23 24 ……… 100]

3. T= 0:1/pi: 1

Therefore T= [0, 0.3183, 0.6366, 0.9549]

4. zeros (1, 3)

The above instruction creates a vector of one row and three columns whose values are
zero

Output= [0 0 0]

5. zeros( 2,4)

Output = 0000

0000

7
9. ones (5,2)

The above instruction creates a vector of five rows and two columns

Output = 11

11

11

11

11

6. a = [ 1 2 3] b = [4 5 6]

a.*b = [4 10 18]

8 if C= [2 2 2]
b.*C results in [8 10 12]

10. plot (t, x)

Ifx = [6 7 8 9] t = [1 2 3 4]

This instruction will display a figure window which indicates the plot of x versus t

8
11. stem (t,x) :- This instruction will display a figure window as shown

12. Subplot: This function divides the figure window into rows and columns.

Subplot (2 2 1) divides the figure window into 2 rows and 2 columns 1 represent number
of the figure

Subplot (3 1 2) divides the figure window into 3 rows and 1 column 2 represent number
of the figure

13. Conv

Syntax: w = conv(u,v)

Description: w = conv(u,v) convolves vectors u and v. Algebraically, convolution is the same

operation as multiplying the polynomials whose coefficients are the elements of u and v.

14. Disp

Syntax: disp(X)

Description: disp(X) displays an array, without printing the array name. If X contains a text
string, the string is displayed.Another way to display an array on the screen is to type its name,
9
but this prints a leading "X=," which is not always desirable.Note that disp does not display
empty arrays.

15. xlabel

Syntax: xlabel('string')

Description: xlabel('string') labels the x-axis of the current axes.

16. ylabel

Syntax : ylabel('string')

Description: ylabel('string') labels the y-axis of the current axes.

17. Title

Syntax : title('string')

Description: title('string') outputs the string at the top and in the center of the current axes.

18. grid on

Syntax : grid on

Description: grid on adds major grid lines to the current axes.

19. FFT Discrete Fourier transform.

FFT(X) is the discrete Fourier transform (DFT) of vector X. For matrices, the FFT operation is
applied to each column. For N-D arrays, the FFT operation operates on the first non-singleton
dimension.
FFT(X,N) is the N-point FFT, padded with zeros if X has less than N points and truncated
if it has more.
20. ABS Absolute value.
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.
21. ANGLE Phase angle.
ANGLE(H) returns the phase angles, in radians, of a matrix with complex elements.

22. INTERP Resample data at a higher rate using lowpass interpolation.

Y = INTERP(X,L) resamples the sequence in vector X at L times the original sample rate.
10
The resulting resampled vector Y is L times longer, LENGTH(Y) = L*LENGTH(X).

23. DECIMATE Resample data at a lower rate after lowpass filtering.

Y = DECIMATE(X,M) resamples the sequence in vector X at 1/M times the original
sample rate. The resulting resampled vector Y is M times shorter, i.e., LENGTH(Y) =
CEIL(LENGTH(X)/M). By default, DECIMATE filters the data with an 8th order Chebyshev
Type I lowpass filter with cutoff frequency .8*(Fs/2)/R, before resampling.

11
Experiment 1: Impulse and step response of the discrete systems

Find out the impulse and step response of the systems described as
i) y(n)  2(0.4)n u(n)  (0.2)n u(n)
ii) y(n)  0.447(1.618) n u(n)  0.447(0.618) n u(n)

The impulse response of a system is the output of the system when the input is an impulse, δ(t), and all
initial conditions are zero. Definition: The step response of a system is the output of the system when the
input is a step, H(t), and all initial conditions are zero.

Program:
i) y(n)  2(0.4)n u(n)  (0.2)n u(n)
clc;
close all;
clear all;
N=40
num=[1 0 0];
den=[1 -0.6 0.08];
[h,t] = impz(num,den,N);
stem(h);

figure;
sys = tf(num,den)
y=step(sys);

stem(y)

12
ii) y(n)  0.447(1.618) n u(n)  0.447(0.618) n u(n)
clc;
close all;
clear all;
N=40
num=[0.2762 -.7236 0];
den=[1 -1 -0.9999];
[h,t] = impz(num,den,N);
stem(h);

figure;
sys = tf(num,den)
y=step(sys);

stem(y)

13
Experiment 2 : VERIFICATION OF SAMPLING THEOREM

Sampling: Is the process of converting a continuous time signal into a discrete time signal. It is the first
step in conversion from analog signal to digital signal.
Sampling theorem: Sampling theorem states that “Exact reconstruction of a continuous time base-band
signal from its samples is possible, if the signal is band-limited and the sampling frequency is greater
than twice the signal bandwidth” i.e. fs > 2W, where W is the signal bandwidth.

Nyquist Rate Sampling: The Nyquist rate is the minimum sampling rate required to avoid aliasing, equal
to the highest modulating frequency contained within the signal. In other words, Nyquist rate is equal to
two sided bandwidth of the signal (Upper and lower sidebands). To avoid aliasing, the sampling rate
must exceed the Nyquist rate. i.e. fs > fN.

Program:

% Sine signal and aliasing

fy=3; %signal frequency in Hz
wy=2*pi*fy; %signal frequency in rad/s
% good sampling frequency
clc;
fs=100; %sampling frequency in Hz
tiv=1/fs; %time interval between samples;
t=0:tiv:(3-tiv); %time intervals set
y=sin(wy*t); %signal data set
subplot(4,1,1); plot(t,y,'k'); %plots figure
axis([0 3 -1.5 1.5]);
title('3Hz sine signal');
ylabel('fs=100');
% too slow sampling frequency
fs=4; %sampling frequency in Hz
tiv=1/fs; %time interval between samples;
t=0:tiv:(3-tiv); %time intervals set
y=sin(wy*t); %signal data set
subplot(4,1,2); plot(t,y,'-kd'); %plots figure
axis([0 3 -1.5 1.5]);
xlabel('seconds');
ylabel('fs=4');

% too slow sampling frequency

fs=10; %sampling frequency in Hz
tiv=1/fs; %time interval between samples;
t=0:tiv:(3-tiv); %time intervals set
y=sin(wy*t); %signal data set
subplot(4,1,3); plot(t,y,'-kd'); %plots figure
axis([0 3 -1.5 1.5]);
xlabel('seconds');
ylabel('fs=10');

14
% too slow sampling frequency
fs=60; %sampling frequency in Hz
tiv=1/fs; %time interval between samples;
t=0:tiv:(3-tiv); %time intervals set
y=sin(wy*t); %signal data set
subplot(4,1,4); plot(t,y,'-kd'); %plots figure
axis([0 3 -1.5 1.5]);
xlabel('seconds');
ylabel('fs=60');

Output:

15
Experiment 3: Generation of signals

Generate different signals using MATLAB

a. Sine wave
b. Cosine wave
c. Unit impulse
d. Unit step wave
e. Square wave
f. Exponential waveform
g. Sawtooth waveform
h. Amplitude-modulated Signal
i. Frequency-modulated signal
Program:
1. Unit impulse signal

clc;
clear all;
close all;
dis p(' UNIT I MP UL S E SI GNAL ') ;
N= input (' Enter Number of Samples: ');
n= - N:1:N
x=[zer os( 1,N) , ze ros( 1,N)]
stem( n, x);
xlabel ('Time') ;
ylabel (' Amplitude'); title('Impulse Signal' );

Output:-
UNIT I MP UL S E S I GNAL
Enter Number of Samples: 6

16
2. Unit step signal

clc;
clear all;
close all;
disp (' UNIT STEP SIGNAL') ;
N=input (' Enter Number of Samples: ');
n= - N:1:N
x= [zer os(1,N) 1 ones (1,N)]
stem( n, x);
xlabel ('Time') ;
ylabel (' Amplitude');
title(' Unit Step Response');

Output:-

UNIT STE P SIGNAL

Number of Samaples : 6

17
3. Unit ramp signal
clc;
clear all;
close all;

disp (' UNIT RAMP SIGNAL') ;

N= input (' Enter Number of Samples : ');
a=input('Enter Amplitude : ')
n= 0:1:N
x=a *n
stem ( n, x) ;
xlabel (' Time') ;
ylabel (' Amplitude');
title (' Unit Ramp Response') ;

Output:-
UNI T RA MP S I GNAL
Enter Number of Samples: 6
Enter Amplitude: 20

18
4. Exponential decaying signal
clc;
clear all; www.jntuworld.com
close all;
disp ('EXPONENTIAL DECAYING SIGNAL') ;
N= input (' Enter Number of Samples: ');
a= 0. 5
n= 0: 0.1:N
x=a. ^ n
stem ( n, x) ;
xlabel ('Time') ;
ylabel ('Amplitude');
title ('Exponential Decaying Signal Response');

Output: -

EXPONE NTIAL DECYI NG SIGNAL

Enter Number of Samples: 6

19
5. Exponential growing signal
clc;
clear all;
close all;
disp ('EXPONE NTIAL GROWING SIGNAL ');

N= input (' E enter Number of Samples: ');

a= 0. 5
n= 0: 0.1:N
x=a. ^ -n
stem (n,x);
xlabel ('Time') ;
ylabel (' Amplitude');
title ('Exponential Growing Signal Response');

Output:-
E XPON NTIAL GROWING SIGNAL
Number of Samaples : 6

20
6. Cosine signal

clc;
clear all;
close all;

disp (' COSI NE SIGNAL ');

N= input (' Enter Number of Samples: ');
n= 0: 0. 1:N
x=cos (n)
stem (n, x);
xlabel ('Time') ;
ylabel ('Amplitude');

title ('Cosine Signal');

Output:-
COSINE SIGNAL
Enter Number of Samaples : 16

21
7. Sine signal

dis p ('SINE SIGNAL ');

N= input (' Enter Number of Samples: ');

n= 0:0. 1:N
x=sin( n) ;
stem (n,x);
xlabel ('Time') ; ylabel (' Amplitude'); title('sine Signal');

Output:-
SINE SIGNAL
Enter Number of Samples : 16

Result:- Thus the MATLAB program for generation of all basic signals was performed and the

22
www.jntuworld.com
output was verified.

Sawtooth Signal generation

%sawtooth signal to be analyzed
clc;
close all;
fy=1; %signal frequency in Hz
wy=2*pi*fy; %signal frequency in rad/s
Ty=1/fy; %signal period in seconds
N=256;
fs=N*fy; %sampling frequency in Hz
tiv=1/fs; %time interval between samples;
t=0:tiv:((3*Ty)-tiv); %time intervals set (3 periods)
y3=sawtooth(wy*t); %signal data set
plot(t,y3,'k');
xlabel('seconds'); title('sawtooth signal (3 periods)');

Output:

23
Generation of amplitude modulated signal

Amplitude modulation is a process by which the wave signal is transmitted by modulating the amplitude
of the signal. It is often called AM and is commonly used in transmitting a piece of information through
a radio carrier wave. Amplitude modulation is mostly used in the form of electronic communication.

Program:

clear all; close all;

clc;
t = 0:0.001:5; %time.
fm = 1;%frequency of message signal.
fc = 10;%frequency of carrier signal.
fs=100*fc;%sampling frequency.
Am = 5;%Amplitude of message signal.
Ac = 5;%Amplitude of carrier signal.
msg =Am.*cos(2*pi*fm*t);%message signal.
carrier = Ac.*cos(2*pi*fc*t);%carrier signal.
%% DSB SC MODULATION AND DEMODULATION.
%===========DSB SC IN TIME DOMAIN==================
dsb_sc = msg.*carrier; %dsbsc modulated wave
%=====DSB SC DEMODULATION TIME DOMAIN============
pmo = 2*dsb_sc.*carrier; %product modulator output
pmo = pmo/Ac;
nf = fm/fs; %normalised frequency
[num, den] = butter(5,3*nf); %butter worth lpf of 5th order
msg_r = filter(num,den,pmo); %demodulated signal after passing through lpf

%================ PLOTTING %%=========================

subplot(4,1,1);
plot(t, msg);
title("MESSAGE SIGNAL (TIME DOMAIN)");
xlabel('time (sec)');
ylabel('amplitude');
gridon;
subplot(4,1,2);
plot(t, carrier);
title("CARRIER SIGNAL (TIME DOMAIN)");
xlabel('time (sec)');
ylabel('amplitude');
gridon;
subplot(4,1,3);
plot(t, dsb_sc);

24
title("MODULATED DSB SC SIGNAL (TIME DOMAIN)");
xlabel('time (sec)');
ylabel('amplitude');
gridon;
subplot(4,1,4);
plot(t, msg_r);
title("DEMODULATED DSB SC SIGNAL (TIME DOMAIN)");
xlabel('time (sec)');
ylabel('amplitude');
gridon;

25
Frequency modulated signal using.
Frequency modulation is a technique or a process of encoding information on a particular signal
(analogue or digital) by varying the carrier wave frequency in accordance with the frequency of the
modulating signal.

Program:

fm=25;
fc=400;
B=10
t=0:0.0001:0.5;
m=cos(2*pi*fm*t);
subplot(3,1,1);
plot(t,m);
xlabel('Time');
ylabel('Amplitude');
title('Message Signal');
grid on;

c=cos(2*pi*fc*t);
subplot(3,1,2);
plot(t,c);
xlabel('Time');
ylabel('Amplitude');
title('Carrier Signal');
grid on;

y=cos(2*pi*fc*t+(B*sin(2*pi*fm*t)));
subplot(3,1,3);
plot(t,y);
xlabel('Time');
ylabel('Amplitude');
title('FM Signal');
grid on;

Output:

26
27
Experiment: 4 Cross-correlation of two signals.

Correlation: Correlation determines the degree of similarity between two signals. If the signals are
identical, then the correlation coefficient is 1; if they are totally different, the correlation coefficient is 0,
and if they are identical except that the phase is shifted by exactly 1800(i.e. mirrored), then the
correlation coefficient is -1.

Autocorrelation: The Autocorrelation of a sequence is correlation of a sequence with itself. The

autocorrelation of a sequence x(n) is defined by,

% Program for computing cross-correlation of the sequences

%x=[1, 2, 3, 4] and h=[4, 3, 2, 1]

clc;
clearall;
closeall;
x=input('enter the 1st sequence');
h=input('enter the 2nd sequence');
y=xcorr(x,h);
figure;subplot(3,1,1);
stem(x);ylabel('Amplitude --.');
xlabel('(a) n --.');
subplot(3,1,2);
stem(h);ylabel('Amplitude --.');
xlabel('(b) n --.');
subplot(3,1,3);
stem(fliplr(y));ylabel('Amplitude --.');
xlabel('(c) n --.');
disp ('The resultant signal is');fliplr(y)

28
Auto Correlation for the given sequence
% Auto Correlation for the given sequence
clc;
close all;
clear all;
x=input('Enter the sequence 1: ');
y=xcorr(x,x);
figure;
subplot(2,1,1);
stem(x);
ylabel('Amplitude->');
xlabel('n->');
title('Input sequence');
subplot(2,1,2);
stem(fliplr(y));
ylabel('amplitude');
xlabel('n->');
title('Output sequence');
disp('the resultant is ');
fliplr(y);

Input: Enter the sequence 1:

[1 2 4 3]

Output:

29
30
Experiment 5: Implementation of FFT

Theory:- A fast Fourier transform (FFT) is an efficient algorithm to compute the discrete Fourier
transform (DFT) and its inverse. FFTs 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 sums of DFT directly would take O(N 2) arithmetical operations. An FFT
is an algorithm to compute the same result in only O(N log N) operations. In general, such
algorithms depend upon the factorization of N, but there are FFTs with O(N log N) complexity
for all N, even for prime N. Since the inverse DFT is the same as the DFT, but with the opposite
sign in the exponent and a 1/N factor, any FFT algorithm can easily be adapted for it as well.

31
Program:
clear all;
N=8;
m=8;
a=input('Enter the input sequence');
n=0:1:N-1;

subplot(2,2,1);
stem(n,a);
xlabel('Time Index n');
ylabel('Amplitude');
title('Sequence');

x=fft(a,m);
k=0:1:N-1;
subplot(2,2,2);
stem(k,abs(x));
ylabel('magnitude');
xlabel('Frequency Index K');
title('Magnitude of the DFT sample');
subplot(2,2,3);
stem(k,angle(x));
xlabel('Frequency Index K');
ylabel('Phase');
title('Phase of DFT sample');
ylabel('Convolution');

Output:-

Enter the input sequence[1 1 1 1 0 0 0 0]

32
www.jntuworld.com
Result:- Thus Fast Fourier Transform is Performed using Matlab.

33
www.jntuworld.com
Discrete Fourier Transform (DFT)

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

The inverse DFT allows us to recover the sequence x[n] from the frequency samples.

X(k) is a complex number (remember ejw=cos(w) + jsin(w). 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‟.

34
clc;
x1 = input('Enter the sequence:');
n = input('Enter the length:');
m = fft(x1,n);
disp('N-point DFT of a given sequence:');
disp(m);
N = 0:1:n-1;

subplot(2,2,1);
stem(N,m);
xlabel('Length');
ylabel('Magnitude of X(k)');
title('Magnitude spectrum:');

an = angle(m);
subplot(2,2,2);
stem(N, an);
xlabel('Length');
ylabel('Phase of X(k)');
title('Phase spectrum:');

Output:-
Enter the sequence:[1 1 0 0]
Enter the length:4
N-point DFT of a given sequence:
Columns 1 through 3

2.0000 1.0000 - 1.0000i 0

Column 4

1.0000 + 1.0000i

35
Result:- Thus Discrete Fourier Transform is Performed using Matlab.

36
www.jntuworld.com
Experiment 6: Program to obtain Linear and circular Convolution

MATLAB code and Output:

%Linear Convolution

clc;

clear all;

x=input('Enter the sequence 1:');

h=input('Enter the sequence 2:');

y=conv(x,h);

subplot(3,1,1);

stem(x);

ylabel('Amplitude->');

xlabel('N->');

title('Input sequence x')

subplot(3,1,2);

stem(h);

ylabel('Amplitude->');

xlabel('N->');

title('Input sequence h')

subplot(3,1,3);

stem(y);

ylabel('Amplitude->');

xlabel('N->');

title(‘linear Convolution');

37
Output:

>> convolution_ok

Enter the first sequence - [1 2 3 4]

Enter the second sequence - [4 3 2 1]

y =

4 11 20 30 20 11 4

38
Experiment 7: To find the frequency response of a given system (Transfer Function/
Difference equation form)

To find the frequency response of the following difference equation y(n) – 5

y(n–1) = x(n) + 4 x(n–1)

Program:-
b = [1, 4]; %Numerator coefficients
a = [1, -5]; %Denominator coefficients
w = -2*pi: pi/256: 2*pi;
[h] = freqz(b, a, w);
subplot(2, 1, 1), plot(w, abs(h));
xlabel('Frequency \omega'), ylabel('Magnitude'); grid
subplot(2, 1, 2), plot(w, angle(h));
xlabel('Frequency \omega'), ylabel('Phase - Radians'); grid

39
Output:-

Result:- Thus the MATLAB program for Frequency Response of Difference Equation was
performed and the output was verified.

40
www.jntuworld.com
Simulation of an M-point Moving Average Filter to filtering out high-frequency components from a
signal composed of a sum of several sinusoidal signals.
Discrete-Time Systems
For the simulation of causal LTI discrete-time systems the command filter can be used. There are
several versions of this command. If we denote
num= [p0 p1 . . . pM],
den= [d0 d1 . . . dN],
then y = filter(num,den,x) generates an output vector y of the same length as the specified input vector x
with zero initial conditions, that is, y[-1] y[-2] = ... =y[-N] = 0. The output can also be computed using y
= filter(num,den,x,ic) where ic = [y[-1], y[-2], ..., y[-N]] is the vector of initial conditions. Access to final
conditions is obtained using [y,fc] filter(num,den,x, ic).

The Moving Average System

A causal version of the three-point smoothing filter is obtained by simply delaying the output by one
sample period, resulting in the FIR filter described by
1
y(n)  ( x[n]  x[n  1]  x[n  2])
3

Generalizing the above equation we obtain

1 M 1
y(n)   x[n  k ]
M k 0
which defines a causalM-point smoothing FIR filter. The system of the above Eq is also known as a
moving average filter . We illustrate its use in filtering high-frequency components from a signal
composed of a sum of several sinusoidal signals.]

% Simulation of an M-point Moving Average Filter

% Generate the input signal
clc;
clearall;
closeall;
n = 0:100;
s1 = cos(2*pi*0.05*n); % A low frequency sinusoid
s2 = cos(2*pi*0.17*n); % A high frequency sinusoid
x = s1+s2;
% Implementation of the moving average filter
M = input('Desired length of the filter = ');
num = ones(1,M);
y = filter(num,1,x)/M;

41
% Display the input and output signals
clf;
subplot(2,2,1);
plot(n,s1);
axis([0, 100, -2, 2]);
xlabel('Time index n'); ylabel('Amplitude');
title('Signal # 1');
subplot(2,2,2);
plot(n,s2);
axis([0, 100, -2, 2]);
xlabel('Time index n'); ylabel('Amplitude');
title('Signal # 2');
subplot(2,2,3);
plot(n,x);
axis([0, 100, -2, 2]);
xlabel('Time index n'); ylabel('Amplitude');
title('Input Signal');
subplot(2,2,4);
plot(n,y);
axis([0, 100, -2, 2]);
xlabel('Time index n'); ylabel('Amplitude');
title('Output Signal');
axis;

Cascade of LTI Systems

In practice a causal LTI discrete-time system of higher order is implemented as a cascade of
lower order causal LTI discrete-time systems. For example, the fourth-order discrete-time
system given below

can be realized as a cascade of two second-order discrete-time systems:

]
6. Generate a sequence x[n], and then uses it as theinput of the fourth-order system, generating the output
y[n]. It then applies the same inputx[n] to Stage No. 1 and finds its output sequence y1[n]. Next, it uses
y1[n] as the inputof Stage No. 2 and finds its output y2[n].

42
clf;
x = [1 zeros(1,40)];% Generate the input
n = 0:40;
% Coefficients of 4th-order system
den = [1 1.6 2.28 1.325 0.68];
num = [0.06 -0.19 0.27 -0.26 0.12];
% Compute the output of 4th-order system
y = filter(num,den,x);
% Coefficients of the two 2nd-order systems
num1 = [0.3 -0.2 0.4];den1 = [1 0.9 0.8];
num2 = [0.2 -0.5 0.3];den2 = [1 0.7 0.85];
% Output y1[n] of the first stage in the cascade
y1 = filter(num1,den1,x);
% Output y2[n] of the second stage in the cascade
y2 = filter(num2,den2,y1);
% Difference between y[n] and y2[n]
d = y - y2;
% Plot output and difference signals
subplot(3,1,1);
stem(n,y);
ylabel('Amplitude');
title('Output of 4th-order Realization');grid;
subplot(3,1,2);
stem(n,y2)
ylabel('Amplitude');
title('Output of Cascade Realization');grid;
subplot(3,1,3);
stem(n,d)
xlabel('Time index n');ylabel('Amplitude');
title('Difference Signal');grid;

43
Write a MATLAB Program to compute the outputs of the above two systems for an Input

clf;
n = 0:299;
x1 = cos(2*pi*10*n/256);
x2 = cos(2*pi*100*n/256);
x = x1+x2;
% Compute the output sequences
num1 = [0.5 0.27 0.77];
y1 = filter(num1,1,x); % Output of System No. 1
den2 = [1 -0.53 0.46];
num2 = [0.45 0.5 0.45];
y2 = filter(num2,den2,x); % Output of System No. 2
% Plot the output sequences
subplot(2,1,1);
plot(n,y1);axis([0 300 -2 2]);
ylabel('Amplitude');
title('Output of System No. 1');grid;
subplot(2,1,2);
plot(n,y2);axis([0 300 -2 2]);
xlabel('Time index n’); ylabel(’Amplitude');
title('Output of System No. 2');grid;

44
45
Experiment 9: Determination of the Power Spectrum of a given signal.

Theory: In statistical signal processing the power spectral density is a positive real function of a
frequency variable associated with a stationary stochastic process, or a deterministic function of
time, which has dimensions of power per Hz, or energy per Hz. It is often called simply the
spectrum of the signal. Intuitively, the spectral density captures the frequency content of a
stochastic process and helps identify periodicities. The PSD is the FT of autocorrelation function,
R(τ) of the signal if the signal can be treated as a wide-sense stationary random process.

Program:-

%Power spectral density

t = 0:0.001:0.6;
x = sin(2*pi*50*t)+sin(2*pi*120*t);
y = x + 2*randn(size(t));
figure,plot(1000*t(1:50),y(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)');

Y = fft(y,512);
%The power spectral density, a measurement of the energy at various frequencies, is:
Pyy = Y.* conj(Y) / 512;
f = 1000*(0:256)/512;
figure,plot(f,Pyy(1:257))
title('Frequency content of y');
xlabel('frequency (Hz)');

46
Output:-

Result:- Thus the MATLAB program for Power Spectral Density was performed and the
output was verified.

47
www.jntuworld.com
Experiment 10: FIR Low pass Filter design
Aim :- To Design FIR LP Filter using Rectangular/Triangular/kaiser Windowing Technique.

Apparatus Used:- S yste m with MAT L A B 6.5.

Algorithm:-
1) Enter the pass band ripple (rp) and stop band ripple (rs).
2) Enter the pass band frequency (fp) and stop band frequency (fs).
3) Get the sampling frequency (f), beta value.
4) Calculate the analog pass band edge frequencies, w1 and
w2. w1 = 2*fp/f
w2 = 2*fs/f
5) calculate the numerator and denominator
6) Use an If condition and ask the user to choose either Rectangular Window or
Triangular window or Kaiser window..
7) use rectwin,triang,kaiser commands
8) Calculate the magnitude of the frequency response in decibels
(dB m=20*log10(abs(h))
9) Plot the magnitude response [magnitude in dB Vs normalized frequency (om/pi)]
10) Give relevant names to x and y axes and give an appropriate title for the plot.
11) Plot all the responses in a single figure window. [Make use of subplot]

48
Program:-

%FIR Filter design window techniques

clc;
clear all;
close all;

rp=input('enter passband ripple');

rs=input('enter the stopband ripple');
fp=input('enter passband freq');
fs=input('enter stopband freq');
f=input('enter sampling freq ');
beta=input('enter beta value„);

wp=2*fp/f;
ws=2*fs/f;
num=-20*log10(sqrt(rp*rs))-13;
dem=14.6*(fs- fp)/f;
n=ceil(num/dem);
n1=n+1;
if(rem(n,2)~=0)
n1=n;
n=n-1;
end

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

if(c==1)
y=rectwin(n1);
disp('Rectangular window filter response');
end
if (c==2)
y=triang(n1);
disp('Triangular window filter response');
end
if(c==3)
y=kaiser(n1,beta);
disp('kaiser window filter response');
end

49
%LPF
b=fir1(n,wp,y);
[h,o]=freqz(b,1,256);
m=20*log10(abs(h));
plot(o/pi,m);
title('LPF');
ylabel('Gain in dB-->');
xlabel('(a) Normalized frequency-->');

Output:-

enter passband ripple 0.02

enter the stopband ripple 0.01
enter passband freq 1000
enter stopband freq 1500
enter sampling freq 10000
enter beta value
enter your choice of window function 1. rectangular 2. triangular 3.kaiser:
1

Rectangular window filter response

50
enter your choice of window function 1. rectangular 2. triangular 3.kaiser:
2
triangular window filter response

enter beta value 5

enter your choice of window function 1. rectangular 2. triangular 3.kaiser:
3
kaiser window filter response

Result:- Thus FIR LP Filter is designed for Rectangular/triangular/kaiser windowing techniques

using MATLAB.

51
Experiment 11: FIR High pass Filter design
To Design FIR HP Filter using Rectangular/Triangular/kaiser Windowing Technique.

Apparatus Used:- S yste m with MAT L A B 6.5.

Algorithm:-
1) Enter the pass band ripple (rp) and stop band ripple (rs).

2) Enter the pass band frequency (fp) and stop band frequency (fs).

3) Get the sampling frequency (f), beta value.

4) Calculate the analog pass band edge frequencies, w1 and w2.

w1 = 2*fp/f

w2 = 2*fs/f
5) calculate the numerator and denominator
6) Use an If condition and ask the user to choose either Rectangular Window or
Triangular window or Kaiser window..

7) use rectwin,triang,kaiser commands

8) Calculate the magnitude of the frequency response in decibels

(dB m=20*log10(abs(h))

9) Plot the magnitude response [magnitude in dB Vs normalized frequency (om/pi)]

10) Give relevant names to x and y axes and give an appropriate title for the plot.

11) Plot all the responses in a single figure window.[Make use of subplot]

52
Program:-

%FIR Filter design window techniques

clc;
clear all;
close all;

rp=input('enter passband ripple');

rs=input('enter the stopband ripple');
fp=input('enter passband freq');
fs=input('enter stopband freq');
f=input('enter sampling freq ');
beta=input('enter beta value');

wp=2*fp/f;
ws=2*fs/f;
num=-20*log10(sqrt(rp*rs))-13;
dem=14.6*(fs- fp)/f;
n=ceil(num/dem);
n1=n+1;
if(rem(n,2)~=0)
n1=n;
n=n-1;
end

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

if(c==1)
y=rectwin(n1);
disp('Rectangular window filter response');
end
if (c==2)
y=triang(n1);
disp('Triangular window filter response');
end
if(c==3)
y=kaiser(n1,beta);
disp('kaiser window filter response');
end

53
%HPF
b=fir1(n,wp,'high',y);
[h,o]=freqz(b,1,256);
m=20*log10(abs(h));
plot(o/pi,m);
title('HPF');
ylabel('Gain in dB-->');
xlabel('(b) Normalized frequency-->');

Output:-

enter passband ripple 0.02

enter the stopband ripple 0.01
enter passband freq 1000
enter stopband freq 1500
enter sampling freq 10000
enter beta value
enter your choice of window function 1. rectangular 2. triangular 3.kaiser:
1

Rectangular window filter response

enter your choice of window function 1. rectangular 2. triangular 3.kaiser:

2
triangular window filter response

54
www.jntuworld.com
enter beta value 5
enter your choice of window function 1. rectangular 2. triangular 3.kaiser:
3
kaiser window filter response

Result:- Thus FIR HP Filter is designed for Rectangular/triangular/kaiser windowing techniques

using MATLAB.

55
www.jntuworld.com
Experiment 12: IIR Low pass Filter design
To Design and generate IIR Butterworth Analog LP Filter using MATLAB

Algorithm:-
1) Enter the pass band ripple (rp) and stop band ripple (rs).

2) Enter the pass band frequency (fp) and stop band frequency (fs).

3) Get the sampling frequency (f).

4) Calculate the analog pass band edge frequencies, w1 and w2.

w1 = 2*fp/f

w2 = 2*fs/f
5) Calculate the order and 3dB cutoff frequency of the analog filter. [Make use of the following
function]
[n,wn]=buttord(w1,w2,rp,rs,‟s‟)

6) Design an nth order analog lowpass Butter worth filter using the following statement.
[b,a]=butter(n,wn,‟s‟)

7) Find the complex frequency response of the filter by using „freqs( )‟

function

[h,om]=freqs(b,a,w) where, w = 0:.01:pi

This function returns complex frequency response vector „h‟ and frequency vector „om‟ in
radians/samples of the filter.

8) Calculate the magnitude of the frequency response in decibels

(dB m=20*log10(abs(h))

9) Plot the magnitude response [magnitude in dB Vs normalized frequency (om/pi)]

10) Calculate the phase response using an = angle(h)
11) Plot the phase response [phase in radians Vs normalized frequency (om/pi)]

12) Give relevant names to x and y axes and give an appropriate title for the plot.

13) Plot all the responses in a single figure window.[Make use of subplot]

56
www.jntuworld.com
Program :-
% IIR filters
clc;
clear all;
close all;
warning off;

disp('enter the IIR filter design specifications');

rp=input('enter the passband ripple');
rs=input('enter the stopband ripple');
wp=input('enter the passband freq');
ws=input('enter the stopband freq');
fs=input('enter the sampling freq');

w1=2*wp/fs;w2=2*ws/fs;
[n,wn]=buttord(w1,w2,rp,rs,'s'); % Find the order n and cutt off frequency

disp('Frequency response of IIR HPF is:');

[b,a]=butter(n,wn,'low','s'); % Find the filter co-efficients of LPF

w=0:.01:pi;
[h,om]=freqs(b,a,w); % Plot the frequency response

m=20*log10(abs(h));
subplot(2,1,1);
plot(om/pi,m);
title('magnitude response of IIR Low Pass filter is:');
xlabel('(a) Normalized freq. -->'); ylabel('Gain in dB-->');

an=angle(h);
subplot(2,1,2);
plot(om/pi,an);
title('phase response of IIR Low Pass filter is:');
xlabel('(b) Normalized freq. -->');
ylabel('Phase in radians-->');

57
www.jntuworld.com
Output:-
enter the IIR filte r design s pecifications
enter the pass band ripple 0. 15
enter the stopband ripple 60
enter the pass band freq1500
enter the stopband freq3000
enter the sampling frq7000
Frequency response of II R L P F is :

Result:- Thus IIR Low Pass Filter is designed using MATLAB.

58
Experiment 13: IIR High pass Filter design
Aim: -To Design and generate IIR Butterworth Analog HP Filter using MATLAB
Algorithm:-
1) Enter the pass band ripple (rp) and stop band ripple (rs).

2) Enter the pass band frequency (fp) and stop band frequency (fs).

3) Get the sampling frequency (f).

4) Calculate the analog pass band edge frequencies, w1 and w2.

w1 = 2*fp/f

w2 = 2*fs/f
5) Calculate the order and 3dB cutoff frequency of the analog filter. [Make use of the following
function]
[n,wn]=buttord(w1,w2,rp,rs,‟s‟)

6) Design an nth order analog high pass Butter worth filter using the following statement.
[b,a]=butter(n,wn,‟high‟,‟s‟)

7) Find the complex frequency response of the filter by using „freqs( )‟ function

[h,om]=freqs(b,a,w) where, w = 0:.01:pi

This function returns complex frequency response vector „h‟ and frequency vector „om‟ in
radians/samples of the filter.

8) Calculate the magnitude of the frequency response in decibels

(dB m=20*log10(abs(h))

9) Plot the magnitude response [magnitude in dB Vs normalized frequency (om/pi)]

10) Calculate the phase response using an = angle(h)

11) Plot the phase response [phase in radians Vs normalized frequency (om/pi)]

12) Give relevant names to x and y axes and give an appropriate title for the plot.

13) Plot all the responses in a single figure window.[Make use of subplot]

59
Program :-

% IIR filters
clc;
clear all;
close all;
warning off;

disp('enter the IIR filter design specifications');

rp=input('enter the passband ripple');
rs=input('enter the stopband ripple');
wp=input('enter the passband freq');
ws=input('enter the stopband freq');
fs=input('enter the sampling freq');

w1=2*wp/fs;w2=2*ws/fs;
[n,wn]=buttord(w1,w2,rp,rs,'s'); % Find the order n and cutt off frequency

disp('Frequency response of IIR HPF is:');

[b,a]=butter(n,wn,'high','s'); % Find the filter co-efficients of HPF

w=0:.01:pi;
[h,om]=freqs(b,a,w); % Plot the frequency response

m=20*log10(abs(h));
subplot(2,1,1);
plot(om/pi,m);
title('magnitude response of IIR High Pass filter is:');
xlabel('(a) Normalized freq. -->'); ylabel('Gain in dB-->');
an=angle(h);
subplot(2,1,2);
plot(om/pi,an);
title('phase response of IIR High Pass filter is:');
xlabel('(b) Normalized freq. -->');
ylabel('Phase in radians-->');

60
Output:-
enter the II R filter design specifications enter the pass
band ripple 0. 15
enter the stop band ripple 60 enter the
pass band freq1500 enter the stop band
freq3000 enter the sampling
freq7000
Frequency response of II R HP F is:

61
Experiment 14: Butterworth filter to square signal
The Butterworth filter is a type of signal processing filter designed to have a frequency response that is
as flat as possible in the pass band. It is also referred to as a maximally flat magnitude filter.
The frequency response of the nth order Butterworth filter is given as

Where ‘n’ indicates the filter order, ‘ω’ = 2πƒ, Epsilon ε is maximum pass band gain, (Amax). If we
define Amax at cut-off frequency -3dB corner point (ƒc), then ε will be equal to one and thus ε2 will also
be equal to one. But, if we want to define Amax at another voltage gain value, consider 1dB, or 1.1220
(1dB = 20logAmax) then the value of ε can be found by:

Where, H0 represents the maximum pass band gain and H1 represents the minimum pass band gain.
Now, if we transpose the above equation, then we will get

Program:

% Response of Butterworth filter to square signal

wc=10; % desired cut-off frequency
N=5; % order of the filter
%analog Butterworth filter:
[num,den]=butter(N,wc,'s');
G=tf(num,den); %transfer function
% Input square signal
fu=1; %signal frequency in Hz
wu=2*pi*fu; %signal frequency in rad/s
fs=100; %sampling frequency in Hz
tiv=1/fs; %time interval between samples;
t=0:tiv:(6-tiv); %time intervals set (6 seconds)
u=square(wu*t); %input signal data set
[y,ty]=lsim(G,u,t); %computes the system output
plot(t,y,'k'); %plots output signal
axis([0 6 -1.5 1.5]);
xlabel('seconds');
title('response to square signal, 5th Butterworth filter');

62
63
Experiment 15: Recovering a sinusoidal signal buried in noise
The most common and obvious problem caused by signal noise is the distortion of the process signal,
causing incorrect interpretation or display of a process condition by the equipment. The addition to
and/or subtraction from the process signal translates into an incorrect process variable.

% Filtering the sine + noise signal

fs =4000; % sampling frequency in Hz
tiv =1/ fs ; % time interval between samples ;
t =0: tiv :(4 - tiv ); % time intervals set (4 seconds )
N= length (t ); % number of data points
yr =0.5 * randn (N ,1); % random signal data set
fy =400; % sinusoidal signal frequency (400 Hz )
ys = sin ( fy *2* pi *t ); % sinusoidal signal
y= ys +yr'; % the signal + noise
% plot sine + noise ( first 0 .1 sec ):
subplot (2 ,1 ,1); plot (t (1:400) , y (1:400) ,'k');
%axis ([0 0 .1 -2 2]);
ylabel ('signal + noise');
title ('filtering the sine + noise signal');
fl =370; % desired low cut - off frequency in Hz
fh =430; % desired high cut - off frequency in Hz
wl = fl *2* pi ; wh = fh *2* pi ; % to rad /s
wb =[ wlwh ]; % the pass band of the filter
N =10; % order of the filter (5+5)
% analog Butterworth filter :
[num , den ]= butter (N ,wb ,'s');
G= tf (num , den ); % transfer function of the filter
yout = lsim (G ,y ,t ); % filter output
% plot extracted signal ( first 0 .1 sec. ):
subplot (2 ,1 ,2); plot (t (1:400) , yout (1:400) ,'k');
%axis ([0 0 .1 -2 2]);
xlabel ('seconds'); ylabel ('extracted signal');
sound (y , fs ); % sound of sine + noise
pause (5);
sound ( yout , fs ); % sound of extracted signal

64
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0
-2
-1
10
signal+noise
2 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1

65
Experiment 16: Implementation of Decimation Process

Theory :-
Sampling rate conversion (SRC) is a process of converting a discrete-time signal at a given
rate to a different rate. This technique is encountered in many application areas such as:
 Digital Audio
 Communications systems
 Speech Processing
 Antenna Systems
 Radar Systems etc
Sampling rates may be changed upward or downward. Increasing the sampling rate is called
interpolation, and decreasing the sampling rate is called decimation. Reducing the sampling rate
by a factor of M is achieved by discarding every M-1 samples, or, equivalently keeping every
M‟th sample. Increasing the sampling rate by a factor of L (interpolation by factor L) is achieved
by inserting L-1 zeros into the output stream after every sample from the input stream of
samples. This system can perform SRC for the following cases:
• Decimation by a factor of M
• Interpolation by a factor ofL
• SRC by a rational factor of L/M.

Decimator :
To reduce the sampling rate by an integer factor M, assume a new sampling period

The re-sampled signal is

The system for performing this operation, called down-sampler, is shown below:

66
Down-sampling generally results in aliasing. Therefore, in order to prevent aliasing, x(n)
should be filtered prior to down-sampling with a low-pass filter that has a cutoff frequency
The cascade of a low-pass filter with a down-sampler illustrated below and is called
decimator.

Program:-

% Illustration of Decimation Process

clc;
close all;
clear all;

M = input('enter Down-sampling factor : ');

N = input('enter number of samples :');
n = 0:N-1;
x = sin(2*pi*0.043*n) + sin(2*pi*0.031*n);
y = decimate(x,M,'fir');

subplot(2,1,1);
stem(n,x(1:N));
title('Input Sequence');
xlabel('Time index n');
ylabel('Amplitude');

subplot(2,1,2);
m = 0:(N/M)-1

www.jntjntuuwororld.ld.ccom
om

67
stem(m,y(1:N /M));
title('Output Sequence');
xlabel('Time index n');ylabel('Amplitude');

Output:-
enter Down-sampling factor : 3
enter number of samples :100

Result:- Thus Decimation Process is implemented using Mat lab.

Note:- Observe the Output Sequence for Different values of M.

68
Experiment 17: Implementation of Interpolation Process
Aim:- To perform interpolation process using Mat lab.
Apparatus required:- S yste m with MAT L A B 6.5.
Theory:- To increase the sampling rate by an integer factor L. If xa(t) is sampled with a
sampling frequency fs = 1/Ts, then

To increase the sampling rate by an integer factor L, it is necessary to extract the samples

from x(n). The samples of xi(n) for values of n that are integer multiples of L are easily extracted
from x(n) as follows:

The system performing the operation is called up-sampler and is shown below:

After up-sampling, it is necessary to remove the frequency scaled images in xi(n), except those
that are at integer multiples of 2π. This is accomplished by filtering xi(n) with a low-pass filter
that has a cutoff frequency of π /L and a gain of L. In the time domain, the low-pass filter inter -
polates between the samples at integer multiples of L as shown below and is called interpolator.

69
Program:-
% Illustration of Interpolation Process
clc;
close all;
clear all;
www.jntjntuuwororld.ld.ccom
om
L = input('Up-sampling factor = ');
N = input('enter number of samples :');
n = 0:N-1;
x = sin(2*pi*0.043*n) + sin(2*pi*0.031*n);
y = interp(x,L);

subplot(2,1,1);
stem(n,x(1:N));
title('Input Sequence');
xlabel('Time index n');
ylabel('Amplitude');

subplot(2,1,2);
m = 0:(N*L)-1;
stem(m,y(1:N*L));
title('Output Sequence');
xlabel('Time index n');
ylabel('Amplitude');

Output:-

70
Result:- Thus Interpolation Process is implemented using Mat lab.

Note:- Observe the Output Sequence for Different values of L.

71
www.jntuworld.com
Experiment- 18. Digital Modulations (ASK,PSK and FSK) of Sine
Signals

Digital Modulation provides more information capacity, high data security, quicker system
availability with great quality communication. Hence, digital modulation techniques have a
greater demand, for their capacity to convey larger amounts of data than analog ones.
There are many types of digital modulation techniques and we can even use a combination of
these techniques as well. In this chapter, we will be discussing the most prominent digital
modulation techniques.

Amplitude Shift Keying

The amplitude of the resultant output depends upon the input data whether it should be a zero
level or a variation of positive and negative, depending upon the carrier frequency.
Amplitude Shift Keying (ASK) is a type of Amplitude Modulation which represents the binary
data in the form of variations in the amplitude of a signal.
Following is the diagram for ASK modulated waveform along with its input.

Any modulated signal has a high frequency carrier. The binary signal when ASK is modulated,
gives a zero value for LOW input and gives the carrier output for HIGH input.

Frequency Shift Keying

The frequency of the output signal will be either high or low, depending upon the input data
applied.
Frequency Shift Keying (FSK) is the digital modulation technique in which the frequency of
the carrier signal varies according to the discrete digital changes. FSK is a scheme of frequency
modulation.
Following is the diagram for FSK modulated waveform along with its input.

72

www.jntuworld.com
The output of a FSK modulated wave is high in frequency for a binary HIGH input and is low in
frequency for a binary LOW input. The binary 1s and 0s are called Mark and Space
frequencies.

Phase Shift Keying

The phase of the output signal gets shifted depending upon the input. These are mainly of two
types, namely BPSK and QPSK, according to the number of phase shifts. The other one is DPSK
which changes the phase according to the previous value.
Phase Shift Keying (PSK) is the digital modulation technique in which the phase of the carrier
signal is changed by varying the sine and cosine inputs at a particular time. PSK technique is
widely used for wireless LANs, bio-metric, contactless operations, along with RFID and
Bluetooth communications.
PSK is of two types, depending upon the phases the signal gets shifted. They are −
Binary Phase Shift Keying (BPSK)
This is also called as 2-phase PSK (or) Phase Reversal Keying. In this technique, the sine wave
carrier takes two phase reversals such as 0° and 180°.
BPSK is basically a DSB-SC (Double Sideband Suppressed Carrier) modulation scheme, for
message being the digital information.
Following is the image of BPSK Modulated output wave along with its input.

73
Program:

ASK:
clc%for clearing the command window
closeall %for closing all the window except command window
clearall %for deleting all the variables from the memory
fc=input('Enter the freq of Sine Wave carrier:'); % ex: 100
fp=input('Enter the freq of Periodic Binary pulse (Message):'); % ex: 10
amp=input('Enter the amplitude (For Carrier & Binary Pulse Message):'); % ex: 4
t=0:0.001:1; % For setting the sampling interval
c=amp.*sin(2*pi*fc*t);% For Generating Carrier Sine wave
subplot(3,1,1) %For Plotting The Carrier wave
plot(t,c)
xlabel('Time')
ylabel('Amplitude')
title('Carrier Wave')
m=amp/2.*square(2*pi*fp*t)+(amp/2);%For Generating Square wave message
subplot(3,1,2) %For Plotting The Square Binary Pulse (Message)
plot(t,m)
xlabel('Time')
ylabel('Amplitude')
title('Binary Message Pulses')
w=c.*m; % The Shift Keyed Wave
subplot(3,1,3) %For Plotting The Amplitude Shift Keyed Wave
plot(t,w)
xlabel('Time')
ylabel('Amplitude')
title('Amplitide Shift Keyed Signal')

FSK:
clc%for clearing the command window
closeall %for closing all the window except command window
clearall %for deleting all the variables from the memory
fc1=input('Enter the freq of 1st Sine Wave carrier:'); % ex: 10

74
fc2=input('Enter the freq of 2nd Sine Wave carrier:'); % ex: 30
fp=input('Enter the freq of Periodic Binary pulse (Message):'); % ex: 5
amp=input('Enter the amplitude (For Both Carrier & Binary Pulse Message):'); % ex: 4
amp=amp/2;
t=0:0.001:1; % For setting the sampling interval
c1=amp.*sin(2*pi*fc1*t);% For Generating 1st Carrier Sine wave
c2=amp.*sin(2*pi*fc2*t);% For Generating 2nd Carrier Sine wave
subplot(4,1,1); %For Plotting The Carrier wave
plot(t,c1)
xlabel('Time')
ylabel('Amplitude')
title('Carrier 1 Wave')
subplot(4,1,2) %For Plotting The Carrier wave
plot(t,c2)
xlabel('Time')
ylabel('Amplitude')
title('Carrier 2 Wave')
m=amp.*square(2*pi*fp*t)+amp;%For Generating Square wave message
subplot(4,1,3) %For Plotting The Square Binary Pulse (Message)
plot(t,m)
xlabel('Time')
ylabel('Amplitude')
title('Binary Message Pulses')
fori=0:1000 %here we are generating the modulated wave
ifm(i+1)==0
mm(i+1)=c2(i+1);
else
mm(i+1)=c1(i+1);
end
end
subplot(4,1,4) %For Plotting The Modulated wave
plot(t,mm)
xlabel('Time')
ylabel('Amplitude')
title('Modulated Wave')

75
PSK:
A=5;
t=0:0.001:1;
f1=20;
f2=5;
x=A.*sin(2*pi*f1*t);
subplot(3,1,1);
plot(t,x);
title('carrier'); grid on;
u=square(2*pi*f2*t);
subplot(3,1,2);
plot(t,u);
title('message signal');
v=x.*u;
subplot(3,1,3);
plot(t,v);
title('PSK'); grid on;

76
77
Experiment: 19. Write a program in MATLAB to compute the
frequency response of the system
(2  z )
x( z ) 
(1  0.6 z)
A spectrum of input sinusoids is applied to a linear time invariant discrete-time system to obtain
the frequency response of the system. The frequency response of the discrete-time system gives
the magnitude and phase response of the system to the input sinusoids at all frequencies.

Functions may be require in the following programs

1.freqz
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
2.impz
impz Impulse response of digital filter
[H,T] = impz (B,A) computes the impulse response of the filter choosing the number of
samples for you, and returns the response in column vector H and a vector of times

clc;

clf;
w = -4*pi:8*pi/511:4*pi;
num = [2 1];den = [1 -0.6];
h = freqz(num, den, w);
% Plot the frequency response
subplot(2,1,1)
plot(w/pi,abs(h));grid
title('Magnitude Spectrum |H(e^{j\omega})|')
xlabel('\omega /\pi');
ylabel('Amplitude');
subplot(2,1,2)
plot(w/pi,angle(h));grid
title('Phase Spectrum arg[H(e^{j\omega})]')
xlabel('\omega /\pi');
ylabel('Phase, radians');

78
Experiment 20: Write a program to find out the impulse
response of any discrete system.
When a system is energizes by a delta function, it produces an output known as its impulse
response. For an LTI system, the impulse response completely determines the output of the
system given any arbitrary input. The output can be found using discrete time convolution.

Program:

clc;
clf;
w = -4*pi:8*pi/511:4*pi;
num = [2 1];den = [1 -0.6];
[h,t] = impz (num,den,w) ;
% Plot the frequency response
subplot(2,1,1)
plot(w/pi,abs(h));grid
title('Magnitude Spectrum |H(e^{j\omega})|')
xlabel('\omega /\pi');
ylabel('Amplitude');
subplot(2,1,2)
plot(w/pi,angle(h));grid
title('Phase Spectrum arg[H(e^{j\omega})]')
xlabel('\omega /\pi');
ylabel('Phase, radians');

79
80
Experiment 21: Design a digital filter which meets the following
specifications using the bilinear method
- sampling frequency: 1 Hz,
- 0 dB attenuation of the DC component,
- maximum attenuation of 1 dB at 0.1 Hz,
- minimum attenuation of 15 dB at 0.15 Hz.
Consider a Butterworth and then a Chebychev model for the filter transfer
function.
% FS: sampling frequency [Hz]
% wp: upper passband edge frequency [rad]
% ws: lower stopband edge frequency [rad]
% atmax: maximum attenuation in the passband[dB]
% atmin: minimum attenuation in the stopband[dB]
% fil: filter type (1 = Butterworth, 0 = Chebyshev)
%[b,a]: numerator/denominator coefficients
Fs = 1; atmax = 1; atmin = 15; wp = 0.1*2*pi;
ws=0.15*2*pi; fil=0; %(or fil=1 according to the filter type)
Ws=2*tan(ws/2);
Wp=2*tan(wp/2);
if fil==1
% Filter design using a Butterworth model
Wnorm=2*log10(Ws/Wp);
N=ceil(log10((10^(0.1*atmin)-1)/(10^(0.1*atmax)-1))/Wnorm);
Whp=Wp/((10^(.1*atmax)-1)^(1/(2*N)));
whp=2*atan(Whp/2);
wn=whp/pi;
[b,a]=butter(N,wn);
else
% Filter design using a Chebyshev1 model
epsilon=sqrt(10^(0.1*atmax)-1);
num=sqrt(10^(0.1*atmin)-1)/epsilon;
N=ceil(acosh(num)/acosh(Ws/Wp))
Whp=Wp*cosh((1/N)*acosh(1/epsilon));
whp=2*atan(Whp/2);
wn=whp/pi;
Rp=atmax;
[b,a]=cheby1(N,Rp,wn);
end
[H,w]=freqz(b,a,512);fq=Fs*w/(2*pi);
subplot(221)
zplane(b,a)
subplot(222)
mag=abs(H);
plot(fq,mag)
ylabel('Amplitude');

81
xlabel('Frequency [Hz]');
axis([0 0.5 0 1.1*max(mag)]); grid
thresh1 = [-atmax*ones(1,length(mag))];
thresh2 = [-atmin*ones(1,length(mag))];
subplot(223)
plot(fq,20*log10(mag),w/pi,thresh1,w/pi,thresh2)
ylabel('Amplitude [dB]');
xlabel('Frequency [Hz]'); grid
axis([0 .5 -1.1*atmin 1]);
xlabel('Frequency [Hz]'); grid
phase_function = unwrap(angle(H));
subplot(224)
plot(fq,phase)
axis([0 .5 1.1*min(phase_function) 1]);
ylabel('Phase [rad]');
xlabel('Frequency [Hz]'); grid

Experiment 22: Design a 25th order FIR lowpass filter with a minimum
transition band, using the window method. Consider a sampling frequency fs
= 20 kHz and the following lower and upper passband frequency edges: f1 = 3
kHz and f2 = 7 kHz. Plot the transfer function, the impulse response and the
zeros of the designed filter.
N=25; f1=3e3; f2=7e3; fs=2e4;
wn1=2*f1/fs; wn2=2*f2/fs;
b=fir1(N,[wn1 wn2],boxcar(N+1));
[h,w] = freqz(b,1); fq = fs*w/(2*pi);
amp=20*log10(abs(h)); phase=unwrap(angle(h))*180/pi;
figure; subplot(221); plot(fq,amp);
axis([0 max(fq) min(amp) 10])
xlabel(‘Frequency [Hz]’); ylabel(‘Amplitude [dB]’) ; grid
subplot(223); plot(fq,phase);
axis([0 max(fq) 1.2*min(phase) 1.2*max(phase)])
xlabel(‘Frequency [Hz]’); ylabel(‘Phase [°]’); grid
subplot(2,2,2); impz(b,1); xlabel(‘n’)
subplot(2,2,4); zplane(b,1);

82
Experiment 23: Design a 24th order FIR highpass filter using the frequency sampling
method. Consider a cutoff frequency fc = 6 kHz and a sampling frequency fs = 20 kHz. Plot
the transfer function, the impulse response and the zeros of the designed filter.

N=24; fc=6e3; fs=20e3; inc=0.01; wcn=2*fc/fs; ff=[0:inc:1];

hh=[zeros(1,wcn/inc-1) ones(1,(1-wcn)/inc+2)];
b=fir2(N,ff,hh); [h,w] = freqz(b,1); fq = fs*w/(2*pi);
amp=20*log10(abs(h)); phase=unwrap(angle(h))*180/pi;
subplot(221); plot(fq,amp); axis([0 max(fq) min(amp) 10])
xlabel(‘Frequency [Hz]’); ylabel(‘Amplitude [dB]’); grid
subplot(223); plot(fq,phase);
axis([0 max(fq) 1.2*min(phase) 1.2*max(phase)])
xlabel(‘Frequency [Hz]’);ylabel(‘Phase [°]’);grid
subplot(2,2,2);impz(b,1);xlabel(‘n’);subplot(2,2,4); zplane(b,1);

Experiment 24: Consider a sum of two sinusoids, whose frequencies are 1 kHz and 1.56
kHz, sampled at 10 kHz. Extract the first sinusoid of this mixture using a lowpass FIR
filter.

fs=10e3; fp=1e3; fc=1.56e3; ondp=.01; ondc=.1; ap=1; ac=0;

[n,fo,mo,w] = remezord([fp fc],[ap ac],[ondpondc],fs);
b = remez(n,fo,mo,w);
f1=1e3;f2=1.56e3; f1n=f1/fs; f2n=f2/fs;
sig1=sin(2*pi*f1n*[0:99]); sig2=sin(2*pi*f2n*[0:99]);
sig=sig1+sig2; y=filter(b,1,sig);
[h,w]=freqz(b,1); fq = fs*w/(2*pi);
amp=20*log10(abs(h));
figure ; plot(fq,amp);
axis([0 max(fq) min(amp) 10])
xlabel(‘Frequency [Hz]’); ylabel(‘Amplitude [dB]’); grid
b.
figure ; subplot(2,2,1);
plot(sig1); title(‘First signal’)
subplot(2,2,3); plot(sig2); title(‘Second signal’)
subplot(2,2,2); plot(sig); title(‘Mixture of the two signals’)
subplot(2,2,4); plot(y); title(‘Filtered signal’)

83