08.508 DSP Lab Manual Part-B

Department of Electronics & Communication, VKCET

08.508 Digital Signal Processing Lab Manual PART-B

Prepared by Assini.H

08.508

Digital Signal Processing Lab Syllabus

Credits: 4 PART A: Experiments on Digital Signal Processors. 1. 2. 3. 4. 5. 6.

Sine wave generation. Real Time FIR Filter implementation (Low-pass, High-pass and Band-pass) Real Time IIR Filter Implementation (Low-pass, High-pass and Band-pass) Pseudo Random Sequence Generator. Real time DFT of sine wave. Sampling a given Analog signal and study of aliasing.

PART B: Experiments on MATLAB or LABVIEW 1. Convolution: Linear Convolution, Circular Convolution, Linear Convolution using Circular Convolution. 2. Random Sequence Generation: Uniform, Rayleigh and Normal Distributions. 3. Discrete Fourier Transform: (Unfolding the spectrum, Frequency Unwrapping). 4. Linear convolution using DFT (Overlap-add and Overlap-Save methods). 5. Design & implementation of IIR filters. (Butterworth and Chebyshev Filters). 6. Design & implementation of FIR filters. (Window method and Frequency sampling Method). 7. Optimal Equiripple Design of FIR filters. 8. Generation of AM, FM & PWM waveforms. 9. Study of sampling rate conversion by a rational factor. 10. Z – Transform.

Department of ECE, VKCET

Page 2

08.508

Digital Signal Processing Lab

Prerequisite knowledge and/or skills: • •

Write computer programs in an object-oriented language such as C or C++. Understand basic linear-system concepts: linear time-shift invariance, system functions, Laplace transforms, z-transforms. Represent digital filters as difference equations, signal flow graphs, z-transform system function, poles and zeroes. Understand Fourier transform properties in continuous and discrete domains. Understand digital filter (IIR and FIR) designing concepts.

• • •

Objectives: • • • • • •

To define and use Discrete Fourier Transforms (DFTs) To design and understand simple finite impulse response (FIR) filters To design and understand simple infinite impulse response (IIR) filters To train the students to design and implement practical DSP systems To program a DSP chip TMS320C6713 to filter signals using Code Composer Studio compiler for the chip. The student should understand how to design algorithms for implementation To introduce signal processing tool boxes in LabVIEW

Outcomes: Students will: a) Have the ability to conduct experiments, as well as to analyze and interpret data in various problems using MATLAB, and DSP starter kit using TMS320C6713 b) Be able to design discrete-time filters and/or implement and verify a filtering system c) Gain knowledge of DSP and be able to design filtering methods in discrete-time domains, as well as to analyze the method in the frequency domain


Page 3

08.508

Digital Signal Processing Lab Index

INTRODUCTION TO MATLAB.............................................................................................................................5 CONVOLUTION USING MATLAB......................................................................................................................14 DFT AND LINEAR CONVOLUTION BY DFT USING MATLAB...............................................................................22 Z-TRANSFORM USING MATLAB.......................................................................................................................36 STUDY OF AM, FM AND PWM WAVEFORM AND RANDOM SEQUENCE

GENERATOR USING MATLAB......45

IIR FILTER DESIGN USING MATLAB..................................................................................................................65 FIR FILTER DESIGN USING MATLAB.................................................................................................................89 OPTIMUM EQUIRIPPLE FIR FILTER USING MATLAB.......................................................................................110 STUDY OF SAMPLING RATE CONVERSION BY A RATIONAL FORM USING MATLAB........................................118 BIBLIOGRAPHY..............................................................................................................................................124


Page 4

08.508 Expt. No.1 Objectives:

Digital Signal Processing Lab INTRODUCTION TO MATLAB a) To familiarize a non user with the fundamentals of MATLAB programming b) To introduce the concept of scripts and functions in MATLAB

Equipments required: PC installed with MATLAB Theory: 1.1 Introduction to MATLAB: What is MATLAB? MATLAB - acronym MATrix LABoratory MATLAB is a software package for computation in engineering, science, and applied mathematics. It offers a powerful programming language, excellent graphics, and a wide range of expert knowledge. MATLAB is published by and a trademark of The Math Works, Inc. Some MATLAB products useful for the lab are: Math and Optimization, Statistics and Data Analysis, Signal Processing and Communications • Signal Processing Toolbox • DSP System Toolbox • Communications System Toolbox • Wavelet Toolbox • Fixed-Point Toolbox • RF Toolbox • Phased Array System Toolbox Application Deployment  MATLAB Compiler How to get into MATLAB? a) Double clicking the Icon

on the Desktop

or b) Using the Start Menu >All Programs>MATLAB The MDE: Matlab Development Environment The MATLAB Environment consists the following a)Command Window b)Workspace and Launch pad tabs c) Command History and Current Directory


Page 5

08.508


1.2 Pre-Requisite For learning MATLAB a) Programming concepts-Procedure oriented b) Matrix Algebra c) Typing Skills and Basic Usage of Computers 1.3 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. Within an M file Programs are of two types a) Matlab script b) Matlab Functions Typing quit in the command prompt>> quit, will close MATLAB MDE 1.4 Data-types in MATLAB There are 15 fundamental data types in MATLAB. Each of these data types is in the form of a matrix or array. This matrix or array is a minimum of 0−by−0 in size and can grow to an n−dimensional array of any size. All of the fundamental data types are shown in lowercase text in the diagram below. Additional data types are user−defined, object−oriented user classes and Java classes. Department of ECE, VKCET

Page 6

08.508


Matlab data types:

1.5 MATLAB Classes There are many different data types, or classes, that you can work with in the MATLAB software. You can build matrices and arrays of floating-point and integer data, characters and strings, and logical true and false states. Matlab classes:

1.6 MATLAB operators: Arithmetic operators: plus - Plus uplus - Unary plus minus - Minus uminus - Unary minus mtimes - Matrix multiply times - Array multiply mpower - Matrix power power - Array power . mldivide - Backslash or left matrix divide Department of ECE, VKCET

+ + * * ^ ^ \ Page 7

08.508 mrdivide ldivide rdivide

Digital Signal Processing Lab - Slash or right matrix divide - Left array divide - Right array divide

/ .\ ./

Relational operators: eq - Equal == ne - Not equal ~= lt - Less than < gt - Greater than > le Less than or equal <= ge - Greater than or equal >= Logical operators: and Logical AND & or - Logical OR | not - Logical NOT ~ xor - Logical EXCLUSIVE OR any True if any element of vector is nonzero all - True if all elements of vector are nonzero 1.7 Getting Started The following commands show how to enter numbers, vectors and matrices, and assign them to variables. >> is the Matlab prompt A sample program- using the plot function -a graph Type in the command Window x=[ 1 2 3 4 5 6 7 8 9 10]; Where the details of x can displayed by typing 'whos x’ in the command window Display Name Size x 1x10

Bytes Class Attributes 80 double

Then type y=sin(x);


Page 8

08.508


Enter plot(x,y);

Instead of typing in the command prompt, repeat the steps in an M file .the M file can be opened from the File menu-. File → New → M file and save with a file name.m Enter the name of the filename.m in the command prompt This is called a script file For any clarification regarding plot etc, which are built in functions type help topic i.e. help plot Entering vectors and matrices; built-in variables and functions; Enter a = 2; a is a constant Enter pi enter x = [1;2;3] x is a column array enter A = [1 2 3;4 5 6;7 8 0] A is a 3x3 matrix A command ‘whos’ will give the details of the all the variables created in the MATLAB workspace


Page 9

08.508


Graphs The simplest graphs to create are plots of points in the Cartesian plane x = [1;2;3;4;5]; y = [0;.25;3;1.5;2]; plot(x,y)

plot(x,y,'r+:') 1 0 .8 0 .6 0 .4 0 .2 0 -0 .2 -0 .4 -0 .6 -0 .8 -1

0

0 .1

0 .2

0 .3

0 .4

0 .5

0 .6

0 .7

0 .8

0 .9

Labeling a Graph To add labels to your graph, the functions xlabel, ylabel, and title can be used as follows: xlabel('x axis') ylabel('Y axis') title('Examples for plotting') E x a m p l e s fo r p lo t t in g 1

0.8 0.6

Y a x is

0.4 0.2

0 -0 . 2

-0 . 4 -0 . 6

-0 . 8

-1

0

0.1

0 .2

0.3

0 .4

0.5

0.6

0 .7

0.8

0.9

x a x is


Page 10

08.508


Example: To plot a composite signal x t =cos 2  100 t  cos 2  200 t  Program: Ts = 0.0002; % sampling time, Ts, is .2 ms n = 1:100; % we will have 100 samples x = cos(2*pi*100*n*Ts) + cos(2*pi*200*n*Ts); plot (n, x, 'b'); xlabel('t'); ylabel('x'); title('x(t) = cos(2 pi 100 t) + cos (2 pi 200 t) '); x (t ) = c o s (2 p i 1 0 0 t ) + c o s (2 p i 2 0 0 t ) 2

1 .5

1

x

0 .5

0

-0 . 5

-1

-1 . 5

0

10

20

30

40

50 t

60

70

80

90

100

1.8 Work to do in the lab: What does the following program do? Programming structures for m=1:10 x(m)=m^2; end; disp(x) If-Else-End Structures a = 5; b = 7; asmall = a < b; % asmall is ``true'' if ~asmall disp('a is greater than b') else disp('a is less than b') end While i = 2; while i<8 i=i+2 end Department of ECE, VKCET

Page 11

08.508


A Word of Caution on Saving Program Files Save any files you create and any MathWorks supplied files that you edit in folders outside of the folder tree in which the MATLAB software is installed. If you keep your files in any of the installed folders, your files can be overwritten when you install a new version of MATLAB. 1.9 Formatting display and graphics commands • format short 3.1416 • format short e 3.1416e+00 • format long 3.14159265358979 • format long e 3.141592653589793e+00 • format compact suppresses extra line feeds (all of the output in this paper is in compact format). • axis( ) sets scaling for the x- and y-axes on the current plot. • grid adds a rectangular grid to the current plot • hold on freezes the current plot so that subsequent graphs will be displayed with the current • hold off releases the current plot; the next plot will erase the current before displaying • subplot puts multiple plots in one graphics window. 1.10 Miscellaneous max(x) returns the largest entry of x, if x is a vector; see help max for the result when x is a k-dimensional array ● min(x) analogous to max ● abs(x) returns an array of the same size as x whose entries are the magnitudes of the entries of x ● size(A) returns a vector with the number of rows, columns, etc. of the kdimensional array A ● length(x) returns the ``length'' of the array, i.e. max(size(A)). ● save fname saves the current variables to the file named fname.mat ● load fname load the variables from the file named fname.mat ● quit exits Matlab 1.11 Examples to do in the lab: 1. Write a program for finding the first 20 Fibonacci numbers. 2. Solve the equations: a) ●

b)

3. Sort a one dimensional array of size N in a) Ascending order b) Descending order and c) find largest and smallest number in the array 4. Simulate dice 5. The voltage across a discharging capacitor is Department of ECE, VKCET

Page 12

08.508


Generate a table and plot voltage v(t) versus time t, for t = 0, 1, 2, ....49 6. Consider complex exponential function Plot the imaginary and real part of the function for t = 0, 0.1, 0.2, ...10 1.12 Function MATLAB functions are similar to C functions or FORTRAN subroutines Function Defintion The first line of a function m-file must be of the following form. function [output_parameter_list] = function_name(input_parameter_list) Example: function sum = addtwo(x,y) % addtwo(x,y) Adds two numbers, vectors, whatever sum = x+y; % print the result = x + y Save this as an M file and execute it from the command prompt e.g. addtwo(5,8) 1.13 Basic Programs to do in the Lab: 1. Write a program to generate the following sequences and plot them using the functions stem and plot: unit sample sequence unit step sequence ramp sequence exponential sequence The input parameter desired by the user is the length of the sequence, delay unit and for exponential sequence A, 2. Write a program to generate square wave and saw tooth wave, plot them using the functions stem and plot. The input data specified by the user are: length of the sequence L, peak value A and period N. For the square wave an additional userspecified parameter is the duty cycle K 3. Write a program to generate a sinusoidal sequence and plot the sequence using the stem and plot functions. The input data specified by the user are desired length L, amplitude A, the angular frequency ω0 and the phase Φ, where and 4. Using the function impz write a program to compute and plot the impulse response of a causal finite discrete time system characterized by difference equation of the form The data to the program are the desired length of impulse response, the constants {pk} and {dk} of the difference equations.


Page 13

08.508


Expt. No.2 Objectives:

CONVOLUTION USING MATLAB a) To study linear convolution without using built in function and the function conv in MATLAB b) To study circular convolution without using built in function in MATLAB

Equipments required: PC installed with MATLAB Theory: Convolution is a formal mathematical operation, just as multiplication, addition, and integration. Addition takes two numbers and produces a third number, while convolution takes two signals and produces a third signal. Convolution is used in the mathematics of many fields, such as probability and statistics. In linear systems, convolution is used to describe the relationship between three signals of interest: the input signal, the impulse response, and the output signal. Convolution is an integral concatenation of two signals. It has many applications in numerous areas of signal processing. The most popular application is the determination of the output signal of a linear time-invariant system by convolving the input signal with the impulse response of the system. The linear convolution of two discrete-time signals x(n) and h(n) is defined by As a system x(n) in input signal, h(n) is impulse response of the signal and y(n) is output signal. Linear convolution obey the properties: a) Commutative property b) Associative property c) Distributive property The circular convolution (or periodic convolution) of two finite signal x(n) and h(n) of length N is defined by

Note: If the two sequences are of different length, the keep the length same by padding zeros. Let x(n) is of length L and h(n) is of length M and M
Page 14

08.508


Explicit circular convolution is rarely used; however implicitly it is used very frequently. Any time DFTs (FFT) or Fourier series are multiplied, there is an underlying circular convolution taking place. The circular convolution result is linear convolution with aliasing error. So by padding zeros to the input sequences by the length of linear convolution sequence aliasing error can be avoid. Let x(n) is of length L and h(n) is of length M, then pad M-1 zeros to the end of the sequence x(n) and L-1 zeros to h(n). Graphical Interpretation of linear convolution: • Reflection of h(k) resulting in h(-k) • Shifting of h(-k) resulting in h(n-k) • Element wise multiplication of the sequences x(k) and h(n-k) • Summation of the product sequence x(k) h(n-k) resulting in the convolution value for y(n) Matrix Interpretation of circular convolution: x(n) as 1xN matrix XN h(n) as NxN matrix HN and in general form

Algorithm: a) Linear convolution without conv: 1. Start 2. Enter th input sequence x(n) and h(n) 3. Compute the length L of x(n) and M of h(n) and find length N = L+M-1 of y(n), where N = L+M-1 4. Pad zero to the end of x(n) and h(n) to avoid the accessing of undefined values 5. Initialize a loop of index i with count value of length N 6. Initialize another loop of index j count value up to the current index of outer loop 7. Perform convolution y(i) = y(i) + x(j) * h(i-j+1) (Matlab doesn't accept zero value indexing, ie h(0)) 8. Continue step 7 up to inner loop count value 9. Continue step 7 and 8 up to the outer loop count value 10. Display the result y(n) 11. Stop


Page 15

08.508


b) Linear convolution with conv: 1. Start 2. Enter the input sequence x(n) and h(n) 3. Use conv function to find linear convolution. (Note: The syntax of conv function is: w = conv(u,v) and description is w = conv(u,v) convolves vectors u and v.) 4. Display the result y(n) 5. Stop c) Circular convolution 1. Start 2. Enter the input sequence x(n) and h(n) 3. Compute the length L of x(n) and M of h(n) and find length N = max(L,M) of y(n) 4. Pad zeros to the end of x(n) and h(n) to keep both sequences as same length 5. Initialize a loop of index i with count value of length N 6. Initialize another loop of index j count value up to N 7. For circular shifting of h(n), compute index for h, k = i-j+1 and if j is negative k =N + k 8. Perform convolution y(i) = y(i) + x(j) * h(k) 9. Continue step 7 up to inner loop count value 10. Continue steps 7 to 9 up to the outer loop count value 11. Display the result y(n) 12. Stop Programs: a) Linear convolution without conv function y = my_conv(x,h) L=length(x); M=length(h); % To pad zeros x1=[x,zeros(1,M-1)]; h1=[h,zeros(1,L-1)]; %Convolution for i=1:(L+M-1) y(i)=0; for j=1:i y(i) =y(i)+x1(j)*h1(i-j+1); end end %plot the inputs and output n=0:length(y)-1; subplot(3,1,1); stem(n,x1); title('LINEAR CONVOLUTION WITHOUT "conv"'); xlabel('n'); ylabel('x(n)') Department of ECE, VKCET

Page 16

08.508

Digital Signal Processing Lab subplot(3,1,2); stem(n,h1); xlabel('n'); ylabel('h(n)') subplot(3,1,3); stem(n,y); xlabel('n'); ylabel('y(n)');

b) Linear convolution with conv: clc; clear all; close all; % input x=input('Enter first sequence in matrix form: '); h=input('Enter second sequence in matrix form: '); %convolution y=conv(x,h); %plotting L=length(x); M=length(h); x=[x,zeros(1,M-1)]; h=[h,zeros(1,L-1)]; n=0:length(y)-1; subplot(3,1,1); stem(n,x); title('LINEAR CONVOLUTION WITH "conv"'); xlabel('n'); ylabel('x(n)') subplot(3,1,2); stem(n,h); xlabel('n'); ylabel('h(n)') subplot(3,1,3); stem(n,y); xlabel('n'); ylabel('y(n)'); c) Circular convolution without built in function: function y = my_circonv(x,h) L=length(x); M=length(h); N=max(L,M); N1=L-M; %zero pading if(N1>=0) h=[h,zeros(1,N1)]; Department of ECE, VKCET

Page 17

08.508

Digital Signal Processing Lab else x=[x,zeros(1,-N1)]; end %circular shifting and convolution for n=1:N y(n)=0; for i=1:N j=n-i+1; if(j<=0) j=N+j; end y(n)=[y(n)+x(i)*h(j)]; end end %plot the inputs and output n=0:N-1; subplot(3,1,1); stem(n,x); xlabel('n'); ylabel('x(n)') subplot(3,1,2); stem(n,h); xlabel('n'); ylabel('h(n)') subplot(3,1,3); stem(n,y); xlabel('n'); ylabel('y(n)');

Results: a) Linear convolution without conv : Inputs: x=[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ]; h=[1 1 1 1 1 1 1 1 1 1 1 1 1 1 ]; Plotted inputs and output: L IN E A R C O N V O L U T IO N W IT H O U T " c o n v"

x (n )

1

0 .5

0

0

5

10

15 n

20

25

30

0

5

10

15 n

20

25

30

0

5

10

15 n

20

25

30

h (n )

1

0 .5

0

y (n )

15 10 5 0


Page 18

08.508


Values displayed: >> x x= 1 1 1 1 1 1 1 1 1 1 >> h h= 1 1 1 1 1 1 1 1 1 1 >> y=my_conv(x,h) y = 1 2 3 4 5 6 7 8 9 10 10 9 8 7 6 5 4 3 2 1

1

1

1

1

1

1

1

1

11

12

13

1

14

1

14

14

13 12

11

b) Linear convolution with conv : Inputs: Enter first sequence in matrix form: [1 2 3 4 5 6 7 6 5 4 3 2 1] Enter second sequence in matrix form: [1 1 1 1 1 1 1 1 1 1 1 1 1] Plotted inputs and output: L IN E A R C O N V O L U T I O N W IT H " c o n v " 8

x (n )

6 4 2 0

0

5

10

15

20

25

15

20

25

15

20

25

n

h (n )

1

0 .5

0

0

5

10 n

60

y (n )

40 20 0

0

5

10 n

Values displayed: >> x x=1 2 3 4 5 6 7 6 5 4 3 2 1 0 0 0 0 0 0 0 0 0 0 0 0 >> h h= 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 >> y y = 1 3 6 10 15 21 28 34 39 43 46 48 49 48 46 43 39 34 28 21 15 10 6 3 1 Department of ECE, VKCET

Page 19

08.508


c) Circular convolution without built in function: Inputs: x= [1 2 3 4 5 6 7 6 5 4 3 2 1];h= [1 1 1 1 1 1 1 1 1 1 1 1 1]; (Same inputs for previous program) Plotted inputs and output: C IR C U L A R C O N V O L U T I O N W IT H O U T B U I L T IN F U N C T IO N 8

x (n )

6 4 2 0

0

2

4

6 n

8

10

12

0

2

4

6 n

8

10

12

0

2

4

6 n

8

10

12

h (n )

1

0.5

0

60

y (n )

40 20 0

Values displayed: >> x x= 1 2 3 4 5 6 7 6 >> h h= 1 1 1 1 1 1 1 1 >> y = my_circonv(x,h) y= 49 49 49 49 49 49 49

5

4

3

2

1

1

1

1

1

1

49

49

49

49

49

49

The output is with aliasing. By linear convolution this can be avoided. Linear convolution by circular convolution: Inputs: x= [1 2 3 4 5 6 7 6 5 4 3 2 1 0 0 0 0 0 0 0 0 0 0 0 0];h= [1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 ];


Page 20

08.508


Plotted output: C IR C U L A R C O N V O L U T IO N W IT H O U T B U IL T IN F U N C T IO N 8

x (n )

6 4 2 0

0

5

10

15

20

25

15

20

25

15

20

25

n

h (n )

1

0 .5

0

0

5

10 n

60

y (n )

40 20 0

0

5

10 n

Values displayed: >>x x=1 2 3 4 5 6 7 6 5 4 3 2 1 0 0 0 0 0 0 0 0 0 0 0 0 >>h h=1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 >>y = my_circonv(x,h) y = 1 3 6 10 15 21 28 34 39 43 46 48 49 48 46 43 39 34 28 21 15 10 6 3 1


Page 21

08.508


Expt. No.3

DFT AND LINEAR CONVOLUTION BY DFT USING MATLAB

Objectives:

a) To study DFT tool without and with fft and ifft in MATLAB b) To perform linear convolution by DFT using MATLAB Equipments required: PC installed with MATLAB 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

Where N is chosen such that L N ≥, where L=length of x[n]. The inverse DFT allows us to recover the sequence x[n] from the frequency samples.

Also the DFT operator

is referred as twiddle factor

X(k) is a complex number (remember ). 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 f s = 8kHz and N=8 point DFT, then in the resulting spectrum, k=1 corresponds to 1kHz frequency. For the same f s 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’. By multiplying two N-point DFTs in the frequency domain, we get the circular convolution in the time domain.

Circular Convolution is made identical to linear convolution by choosing the number of points in DFT as N ≥ length of x + length of h -1. For Circular Convolution, N = max (length of x, length of h). Department of ECE, VKCET

Page 22

08.508


The short signals provide useful illustrations of convolution, but they are not typical of practical filtering problems. In many applications the input x[n] is very long with respect to the filter impulse response h[n]. It is often impractical to process the entire input signal at once. An alternative solution is to break the input into blocks, filter each block separately and then combine the blocks appropriately to obtain the output. There are two basic strategies for block processing, which are known as overlap-add and overlap-save, respectively. Overlap-Add Implementation The overlap-add method works by breaking the long input signal into small nonoverlapping sections. If the length of these sections is L and the length of the impulse response is P, then an N= (L+P-1) - point DFT is used to implement the filter for each section. (The (L+P-1)-point DFT avoids problems with time-aliasing inherent in the circular convolution.) The result for each block overlaps with the following block. These overlapping sections must be added together to compute the output. Overlap-Save Implementation The overlap-save method uses a different strategy to break up the input signal. If the length of the circular convolution (i.e., the DFT length) is chosen to be N=L+P-1 , overlapsave uses input segments of length N. The starting location of each input segment is skipped by an amount L, so there is an overlap of P - 1 point with the previous section. Thus this method could also be called the overlapped inputs method. The L “good” points resulting from the N-point circular convolution of each section with the filter are retained. No additions are needed to create the output. Algorithms: a) DFT without using fft 1. Start 2. Enter the input sequence xn and number of points N 3. Compute the length of xn, L and if L>N discard last L-N points of xn otherwise pad N-L zeros to x(n) 4. Create indices n=0,1,....N-1 and k=0,1,...N-1 kn 5. Create matrix for twiddle factor WN and compute matrix W N kn 6. Compute DFT by matrix multiplication Xk = xn * W N 7. Stop b) Plot the spectrum 1. Start 2. Enter the signal Xk 3. Compute the length of Xk , N and create the index k=0,1,... N-1 4. Compute the |Xk| (Note: Use MATLAB function abs( ). To know this function type help abs in command window and read the features of abs or use MATLAB help) 5. Compute the ∟Xk (Note: MATLAB function angle( ) .To know this function type help angle in command window and read the features of angle or use MATLAB help) Department of ECE, VKCET

Page 23

08.508

Digital Signal Processing Lab 6. Plot |Xk| versus k along with label in x and y axis and title 'Magnitude spectrum of X(k)' 7. Plot ∟Xk versus k along with label in x and y axis and title 'Phase spectrum of X(k)' 8. Stop c) IDFT without using ifft 1. Start 2. Enter the input sequence Xk and number of points N 3. Compute the length of Xk , L and if L>N discard last L-N points of xn otherwise pad N-L zeros to x(n) 4. Create indices n=0,1,....N-1 and k=0,1,...N-1 −kn 5. Create matrix for twiddle factor WN and compute matrix W N 1

−kn

6. Compute IDFT by matrix multiplication xn = N ( Xk * W N ) 7. Stop (Note: After completing programs a and c verify your data with DFT pair) d) DFT using fft Note: The syntax of fft is: fft(x,N) is the N-point FFT, padded with zeros if length of x is less than N points and truncated if it has more. 1. Start 2. Enter the input sequence xn and number of points N 3. Compute DFT, Xk = fft(xn,N) 4. Compute length of xn, L 5. Create index n=0,1,... L-1 and k=0,1,....N-1 6. Display xn and Xk 7. Plot xn versus n, real part of Xk versus k and imaginary part of X k versus k (Note: For real part use real( )and for imaginary part imag( ) MATLAB function) 8. Stop e) IDFT using ifft Note: The syntax of fft is: ifft(X,N) is the N-point FFT, padded with zeros if length of X is less than N points and truncated if it has more. 1. Start 2. Enter the input sequence Xk and number of points N 3. Compute IDFT, xn = ifft(Xk, N) 4. Create indices n=0,1,....N-1 and k=0,1,...N-1 5. Display xn 6. Plot xn versus n, real part of Xk versus k and imaginary part of Xk versus k 7. Stop f) Linear convolution using DFT by overlap-add method 1. Start 2. Enter input xn , hn and length of block L Department of ECE, VKCET

Page 24

08.508

Digital Signal Processing Lab 3. 4. 5. 6. 7.

Break the input signal xn into non-overlapping blocks xmn of length L. Compute length of hn, M Zero pad hn to be of length N = L + M-1 Take N-DFT of hn to give Hk, k = 0,1,....,N-1 For each block m: Zero pad xmn to be of length N = L + M-1 Take N-DFT of xmn to give Xmk, k = 0,1,.....,N-1 Multiply: Ymk = Xmk..* Hk, k = 0,1,.....,N-1 Take N-IDFT of Ymk to give ymn, n = 0,1,......N-1 8. Form yn by overlapping the last M-1 samples of y mn with the first M-1 samples of y(m+1)n and adding the result. 9. Stop g) Linear convolution using DFT by overlap-save method 1. Start 2. Enter input xn , hn and length of block L 3. Compute length of hn, M 4. Insert M-1 zeros at the beginning of the input sequence xn 5. Break the padded input signal into overlapping blocks x mn of length N = L + M 1. 6. Zero pad hn to be of length N = L + M-1 7. Take N-DFT of hn to give Hk, k = 0,1,........,N-1 8. For each block m: Take N-DFT of xmn to give Xmk, k =0,1,........,N-1 Multiply: Ymk = Xmk .* Hk, k = 0,1,........,N-1 Take N-IDFT of Ymk to give ymn, n = 0,1,........,N-1 Discard the first M-1 points of each output block ymn 9. Form yn by appending the remaining (i.e., last) L samples of each block ymn 10. Stop Programs: a) DFT without fft %DFT without fft function [Xk]=dft(x,N) %padding zeros or truncating x(n) points L=length(x); if N>L xn=[x,zeros(1,N-L)]; else for i=1:N xn(i)=x(i); end end %Twiddle factor n=0:1:N-1; k=0:1:N-1; nk=n'*k; %Transpose of n or k is must Department of ECE, VKCET

Page 25

08.508

Digital Signal Processing Lab WN=exp(-j*2*pi/N); WNnk=WN.^nk; %DFT Xk=xn*WNnk;

b) To plot Spectrum %Program to plot spectrum function [ ]=my_spectrum(Xk) % length of Xk N=length(Xk); k=0:1:N; %Magnitude spectrum Xk(N+1)=Xk(1); %Periodic property to get N+1 points MXk=abs(Xk); display(‘Press any key to plot the MAGNITUDE & PHASE SPECTRUM of Xk’); pause; subplot(2,1,1); stem(k,MXk); grid on; title('Magnitude Spectrum of X(k)'); xlabel('k'); ylabel('Mag. of X(k)'); %Phase Spectrum PXk=angle(Xk); subplot(2,1,2); grid on; stem(k,PXk); title('Phase Spectrum of X(k)'); xlabel('k'); ylabel('Phase of X(k)'); c) IDFT without ifft %DFT without fft function [xn]=idft(X,N) %padding zeros or truncating x(n) points L=length(X); if N>L Xk=[X,zeros(1,N-L)]; else for i=1:N Xk(i)=X(i); end end %Twiddle factor n=0:1:N-1; Department of ECE, VKCET

Page 26

08.508

Digital Signal Processing Lab k=0:1:N-1; nk=n'*k; %Transpose of n or k is must WN=exp(-j*2*pi/N); WNnk=WN.^-nk; %DFT xn=(Xk*WNnk)/N;

d) DFT using fft %DFT using 'fft' clc; clear all; close all; %input xn=input('Enter the input in matrix form: '); N=input('Enter the number of points: '); %DFT Xk=fft(xn,N); %Display disp('Press any key to display x(n) and X(k)'); pause; disp('x(n) = '); disp(xn); disp('X(k) '); disp(Xk); %Plot disp('Press any key to plot input and output'); pause; n=0:1:N; k=0:1:N; xn(N+1)=0; Xk(N+1)=Xk(1); subplot(3,1,1); stem(n,xn); title('Input x(n)'); xlabel('n'); ylabel('x(n)'); RXk=real(Xk); %real part of X(k) IXk=imag(Xk); %imaginary part of X(k) subplot(3,1,2); stem(k,RXk); title('X(k)'); xlabel('k'); ylabel('Real part of X(k)'); subplot(3,1,3); stem(k,IXk); title('X(k)'); Department of ECE, VKCET

Page 27

08.508


xlabel('k'); ylabel('Imaginary part of X(k)'); e) IDFT using ifft %DFT using 'fft' clc; clear all; close all; %input Xk=input('Enter the input in matrix form: '); N=input('Enter the number of points: '); %IDFT xn=ifft(Xk,N); %Display disp('Press any key to display X(k) and x(n)'); pause; disp('X(k) = '); disp(Xk); disp('x(n)'); disp(xn); f) Linear convolution using DFT by overlap-add method %Linear convolution by DFT using Overlap-add method function [y] = ovradd(x,h,N) M = length(h); Nx = length(x); if(N
Page 28

08.508

Digital Signal Processing Lab lengths % Overlap-add Method %Blocks of x(n) padded with zeros in the last M-1 positions xadd=zeros(B,N); for i=1:B xadd(i,:)=[xsplit(i,:) zeros(1,M-1)]; end %Circular convolution operation for m=1:B ym(m,:)=cconv(xadd(m,:),hz,N); end %Adding the first M-1 terms of next block to the last M-1 terms of previous block y=zeros(B, N); if(B>1) for i=1:B-1 temp=ym(i+1,:); yi(i,:)=[ym(i,1:end-M+1) (ym(i,end-M+2:end)+temp(1,1:M-1))]; end yi(B,:)=[ym(B,:)]; else yi=ym; end %Output to one-dimensional array y=zeros(1,Nx+M-1); for i=1:B if(i==1) y=yi(i,:); if(B==1) y=[yi(i,1:end-M+1)]; end elseif(i==B) y=[y yi(i,M:end-M+3)]; else y=[y yi(i,M:end)]; end end len_y=len_x+M-1; y=[y(1:len_y)]; disp(y');


Page 29

08.508


g) Linear convolution using DFT by overlap-save method %Linear convolution using DFT by Overlap-save method function [y] = ovrsave(x,h,N) M = length(h); Nx = length(x); len_x=Nx; if(N
Page 30

08.508

Digital Signal Processing Lab %Output to proper form y=y'; y=[y(1:(len_x+M-1))];

Results: a) DFT without fft Input >> x=[1 2 3 4 3 2 1]; >> N=8; Output >> Xk=dft(x,N) Xk = 16.0000 -4.8284 - 4.8284i 0 0.8284 - 0.8284i 0.8284 + 0.8284i 0 -4.8284 + 4.8284i b) To plot Spectrum Input to DFT >> x=[1 1 1 1]; N=4; >> Xk=dft(x,N); To plot spectrum >> my_spectrum(Xk) Press any key to plot the MAGNITUDE & PHASE SPECTRUM of Xk M a g n it u d e S p e c t ru m o f X (k ) M a g . o f X (k )

4 3 2 1 0

0

0 .5

1

0

0 .5

1

1 .5

2 2 .5 k P h a s e S p e c tru m o f X (k )

3

3 .5

4

3

3 .5

4

P h a s e o f X (k )

0

-1

-2

-3

1 .5

2 k

2 .5

>> x=[1 1 1 1]; N=16; >> Xk=dft(x,N); >> my_spectrum(Xk)


Page 31

08.508

Digital Signal Processing Lab Press any key to plot the MAGNITUDE & PHASE SPECTRUM of Xk M a g n it u d e S p e c t ru m o f X (k ) M a g . o f X (k )

4 3 2 1 0

0

2

4

0

2

4

6

8 10 k P h a s e S p e c t ru m o f X (k )

12

14

16

12

14

16


4

2

0

-2

6

8 k

10

x=[1 1 1 1]; N=128; >> Xk=dft(x,N); >> my_spectrum(Xk) Press any key to plot the MAGNITUDE & PHASE SPECTRUM of Xk M a g n it u d e S p e c t r u m o f X ( k )

M a g . o f X (k )

4 3 2 1 0

0

20

40

0

20

40

60

80 k P h a s e S p e c t ru m o f X (k )

100

120

140

100

120

140


4 2 0 -2 -4

60

80 k

c) IDFT without ifft Input >> Xk=[16 -4.8284-4.8284i 0 0.8284-0.8284i 0 0.8284+0.8284i 0-8284+4.8284i]; >> N=8; Output >> xn=idft(Xk,N) xn =1.0000 2.0000 3.0000 - 0.0000i 4.0000 - 0.0000i 3.0000 - 0.0000i 2.0000 1.0000 + 0.0000i 0.0000 + 0.0000i Note: Imaginary part of x(n) is not exact zero due to truncating of Twiddle factor. To avoid this use abs( ) Department of ECE, VKCET

Page 32

08.508

Digital Signal Processing Lab >> xn=abs(idft(Xk,N)) xn = 1.0000 2.0000 3.0000

4.0000

3.0000

2.0000

1.0000

0.0000

d) DFT using fft Inputs Enter the input in matrix form: [1 1 1 1] Enter the number of points: 4 Outputs Press any key to display x(n) and X(k) x(n) = 1 1 1 1 X(k) 4 0 0 0 Press any key to plot input and output In p u t x ( n ) x (n )

1 0 .5

I m a g i n a r y p a r t o f X ( k R) e a l p a r t o f X ( k )

0

0

0 .5

1

1 .5

2 n X (k )

2 .5

3

3 .5

4

0

0 .5

1

1 .5

2 k X (k )

2 .5

3

3 .5

4

0

0 .5

1

1 .5

2 k

2 .5

3

3 .5

4

4 2 0

1 0 -1

Enter the input in matrix form: [1 1 1 1] Enter the number of points: 16 In p u t x ( n ) x (n )

1 0 .5

I m a g in a ry p a rt o f X (k ) R e a l p a rt o f X (k )

0

0

2

4

6

8 n X (k )

10

12

14

16

0

2

4

6

8 k X (k )

10

12

14

16

0

2

4

6

8 k

10

12

14

16

5 0 -5

5 0 -5


Page 33

08.508


Enter the input in matrix form: [1 1 1 1] Enter the number of points: 128 In p u t x ( n )

x (n )

1 0 .5 0

0

20

40

60

80

100

120

140

80

100

120

140

80

100

120

140

I m a g in a ry p a rt o f X (k )

R e a l p a rt o f X (k )

n X (k ) 5 0 -5

0

20

40

60 k X (k )

5 0 -5

0

20

40

60 k

e) IDFT using ifft Input Enter the input in matrix form: [16 -4.8284-4.8284i 0 0.8284-0.8284i 0 0.8284+0.8284i 0 -4.8284+4.8284i] Enter the number of points: 8 Output Press any key to display X(k) and x(n) X(k) 16.0000 -4.8284 - 4.8284i 0 0.8284 - 0.8284i 0.8284 + 0.8284i 0 -4.8284 + 4.8284i x(n) 1.0000 2.0000 3.0000 4.0000 3.0000 2.0000 1.0000 0.0000 f) Linear convolution using DFT by overlap-add method Input >> x=[0:1:15 14:-1:0] x = 0 1 2 3 4 5 6 7 8 9 10 11 12 12 11 10 9 8 7 6 5 4 3 2 1 0 >> h=[1 1 1 1]; >> N=10;


13

14

15

14 13

Page 34

08.508

Digital Signal Processing Lab >> ovradd(x,h,N); 0 1.0000 3.0000 6.0000 10.0000 14.0000 18.0000 22.0000 26.0000 30.0000 34.0000 38.0000 42.0000 46.0000 50.0000 54.0000 56.0000 56.0000 54.0000 50.0000 46.0000 42.0000 38.0000 34.0000 30.0000 26.0000 22.0000 18.0000 14.0000 10.0000 6.0000 3.0000 1.0000 0

g) Linear convolution using DFT by overlap-save method >> x=[0:1:15 14:-1:0] x =0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1 0 >> h=[1 1 1 1]; >> N=20; >> ovrsave(x,h,N); 0 1.0000 3.0000 6.0000 10.0000 14.0000 18.0000 22.0000 26.0000 30.0000 34.0000 38.0000 42.0000 46.0000 50.0000 54.0000 56.0000 56.0000 54.0000 50.0000 46.0000 42.0000 38.0000 34.0000 30.0000 26.0000 22.0000 18.0000 14.0000 10.0000 6.0000 3.0000 1.0000 0


Page 35

08.508 Expt. No.4

Digital Signal Processing Lab Z-TRANSFORM USING MATLAB

Objectives: To study z-transform and inverse z-transform for representing discrete-time systems using MATLAB Equipments required: PC installed with MATLAB Theory: The Z-transform, like many integral transforms, can be defined as either a one-sided or two-sided transform. Bilateral Z-transform The bilateral or two-sided Z-transform of a discrete-time signal x[n] is the function X(z) defined as

where n is an integer, j is the imaginary unit, and z is, in general, a complex number: where A is the magnitude of z, and φ is the complex argument (also referred to as angle or phase) in radians. The set of z values for which X(z) exists is called the region of convergence (ROC) and is given by

For some positive numbers Rx- and Rx+. Unilateral Z-transform Alternatively, in cases where x[n] is defined only for n ≥ 0, the single-sided or unilateral Z-transform is defined as

In signal processing, this definition can be used to evaluate the Z-transform of the unit impulse response of a discrete-time causal system. An important example of the unilateral Z-transform is the probability-generating function, where the component x[n] is the probability that a discrete random variable takes the value n, and the function X(z) is usually written as X(s), in terms of s = z − 1. The properties of Z-transforms have useful interpretations in the context of probability theory.


Page 36

08.508

Digital Signal Processing Lab Properties of the ROC: 1. The ROC is always bounded by a circle since the convergence condition is on the magnitude |z|. 2. The ROC for right-sided sequences is always outside of a circle of radius Rx-. 3. The ROC for left-sided sequences is always inside of a circle of radius Rx+. 4. The ROC for two-sided sequences is always an open ring if it exists. 6. 7. 8. 9.

5. 5. The sequences that are zero for n < n l and n > n2 are called finite-duration sequences. The ROC for such sequences is the entire z-plane. If n l < 0, then z = ∞ is not in the ROC. If n2 > 0, then z = 0 is not in the ROC. The ROC cannot include a pole since X (z) converges uniformly in there. There is at least one pole on the boundary of a ROC of a rational X (z). The ROC is one continuous region; that is, the ROC does not come in pieces.

Consider some examples for ROC and the pole-zero plot of a sequence 1.

2.

ROCs are not specified, hence three possible ROCs are:


Page 37

08.508


Inverse Z-transform The inverse Z-transform is

Where

is a counterclockwise closed path encircling the origin and entirely in the region

of convergence (ROC). The contour or path,

, must encircle all of the poles of

.

The inverse z-transform computation requires an evaluation of a complex contour integral that, in general, is a complicated procedure. The most practical approach is to use the partial fraction expansion method In general

It can be expressed as

Where the fist term on the right-hand side is the proper rational part and the second term is the polynomial (finite length) part. This can be obtained by performing polynomial division if M≥N using the deconv function. To perform a partial fraction expansion on the proper rational part of X(z) to obtain

Where pk is the kth pole of X(z) and Rk is the residue at pk. It is assumed that the poles are distinct for which the residues are given by

For repeated poles for the expansion for X (z) has a more general form. If a pole pk has multiplicity r, then its expansion is given by

Where the residues Rk,l are computed using a more general formula Then


Page 38

08.508


Algorithms: a) z-transform using MATLAB symbolic toolbox 1. Start 2. Define symbols n and z 3. Enter input 4. Compute z-transform by ztrans( ) 5. Display the result 6. Stop Note: the syntax of ztrans is F = ZTRANS(f) is the Z-transform of the scalar sym f with default independent variable n. The default return is a function of z F = ZTRANS(f,w) makes F a function of the sym w instead of the default z F = ZTRANS(f,k,w) takes f to be a function of the sym variable k b) Inverse z-transform using MATLAB symbolic toolbox 1. 2. 3. 4. 5. 6.

Start Define symbols n and z Enter input Compute z-transform by iztrans( ) Display the result Stop

Note: the syntax of iztrans is i) f = IZTRANS(F) is the inverse Z-transform of the scalar sym F with default independent variable z. The default return is a function of n ii) f = IZTRANS(F,k) makes f a function of k instead of the default n. Here k is a scalar sym iii) f = IZTRANS(F,w,k) takes F to be a function of w instead of the default symvar(F) and returns a function of k c) Power series expansion using MATLAB 1. 2. 3. 4. 5.

Start Enter numerator and denominator coefficients in matrix form Compute inverse z-transform using deconv( ) Display the result Stop

Note: The syntax and description for deconv( ) are: 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 If A and B are vectors of polynomial coefficients, deconvolution is equivalent to polynomial division Department of ECE, VKCET

Page 39

08.508

Digital Signal Processing Lab d) Partial fraction expansion using MATLAB 1. 2. 3. 4. 5.

Start Enter numerator and denominator coefficients in matrix form Compute residues, poles and direct form expansion using residuez( ) Display the result Stop

Note: The syntax and description for residuez( )are: RESIDUEZ Z-transform partial-fraction expansion [R,P,K] = RESIDUEZ(B,A) finds the residues, poles and direct terms of the partialfraction expansion of B and A are the numerator and denominator polynomial coefficients, respectively, in ascending powers of . R and P are column vectors containing the residues and poles, respectively. K contains the direct terms in a row vector. The number of poles is n = length(A)-1 = length(R) = length(P) The direct term coefficient vector is empty if length(B) < length(A); otherwise, length(K) = length(B)-length(A)+1 e) Pole-zero plotting from system function (z-domain) 1. 2. 3. 4. 5.

Start Enter numerator and denominator in matrix form Compute zeros and poles using roots( ) Plot zeros and poles using zplane( ) Stop

Note i) The syntax and details of roots ( ) are: ROOTS Find polynomial roots ROOTS(C) computes the roots of the polynomial whose coefficients are the elements of the vector C ii) The syntax and details of zplane ( ) are: ZPLANE Z-plane zero-pole plot ZPLANE(Z,P) plots the zeros Z and poles P (in column vectors) with the unit circle for reference. Each zero is represented with a 'o' and each pole with a 'x' on the plot. Multiple zeros and poles are indicated by the multiplicity number shown to the upper right of the zero or pole. ZPLANE(B,A) where B and A are row vectors containing transfer function polynomial coefficients plots the poles and zeros of B(z)/A(z). Note that if B and A are both scalars they will be interpreted as Z and P.


Page 40

08.508


Programs: a) z-transform using MATLAB symbolic toolbox %z-transform using symbolic math toolbox and ztrans( ) clc; clear all; close all; %Define symbols syms n z w k; %Input xn=input('Enter input x(n) in symbolic form: '); %z-transform Xz=ztrans(xn); Xw=ztrans(xn,w); Xkw=ztrans(xn,k,w); %Display disp('x(n)'); disp(xn); disp('X(z)'); disp(Xz); disp('X(w)'); disp(Xw); disp('Xk(w)'); disp(Xkw); b) Inverse z-transform using MATLAB symbolic toolbox %Inverse z-transform using symbolic math toolbox and iztrans( ) clc; clear all; close all; %Define symbols syms n z w k; %Input Xz=input('Enter input X(z)) in symbolic form: '); %z-transform xn=iztrans(Xz); xw=iztrans(Xz,w); %Display disp('X(z)'); disp(xn); disp('x(n)'); disp(xn); disp('x(w)'); disp(xw);


Page 41

08.508


c) Power series expansion using MATLAB % inverse z-transform using the power series method clc; clear all; close all; % Input num=input('Enter the numerator coefficients of X(z) in matrix form: '); den=input('Enter the denominator coefficients of X(z) in matrix form: '); L=input('Enter the length of the power series expansion: '); %Padding zeros to numerator to keep required length of power series num=[num, zeros(1,L-1)]; % Deconvolution [x,r]=deconv(num,den); % Display disp('x(n)'); disp(x); disp('reminder'); disp(r); d) Partial fraction expansion using MATLAB % inverse z-transform using partial fraction clc; clear all; close all; % Input num=input('Enter the numerator coefficients of X(z) in matrix form: '); den=input('Enter the denominator coefficients of X(z) in matrix form: '); %z-transform Partial fraction expansion [r,p,k]=residuez(num,den); % Display disp('Residues'); disp(r'); disp('Poles'); disp(p'); disp('Direct form'); disp(k'); e) Pole-zero plotting from system function (z-domain) % pole zero plot of z-transform using zplane clc; clear all; close all; % Input num=input('Enter the numerator coefficients of X(z) in matrix form: '); den=input('Enter the denominator coefficients of X(z) in matrix form: '); %Poles and zeros z=roots(num); Department of ECE, VKCET

Page 42

08.508

Digital Signal Processing Lab p=roots(den); % Display zplane(z,p) grid on; title('Pole-Zero plot of X(z)');

Results: a) z-transform using MATLAB symbolic toolbox Given that Enter input x(n) in symbolic form: 1+n+n^2 x(n) 1+n+n^2 X(z) z/(z-1)+z/(z-1)^2+z*(z+1)/(z-1)^3 X(w) w/(w-1)+w/(w-1)^2+w*(w+1)/(w-1)^3 Xk(w) (1+n+n^2)*w/(w-1) b) Inverse z-transform using MATLAB symbolic toolbox Given that Enter input X(z)) in symbolic form: 1/(1-z^-2) X(z) 1/2+1/2*(-1)^n x(n) 1/2+1/2*(-1)^n x(w) 1/2+1/2*(-1)^w c) Power series expansion using MATLAB Given that Enter the numerator coefficients of X(z) in matrix form: [1 0.5 1] Enter the denominator coefficients of X(z) in matrix form: [1 0 1 2] Enter the length of the power series expansion: 5 x(n) 1.0000 0.5000 0 -2.5000 -1.0000 Reminder 0 0 0 0 0 2.5000 6.0000 2.0000 d) Partial fraction expansion using MATLAB 1) Given that Enter the numerator coefficients of X(z) in matrix form: [1 -1/4] Enter the denominator coefficients of X(z) in matrix form: [1 -1/6 -1/6] Residues 0.3000 0.7000 Poles Department of ECE, VKCET

Page 43

08.508

Digital Signal Processing Lab 0.5000 -0.3333 Direct form (Null) 2) Given that

Enter the numerator coefficients of X(z) in matrix form: [1 0 1 .5] Enter the denominator coefficients of X(z) in matrix form: [1 -1/6 0 -1/6] Residues 2.3843 0.8078 - 0.9877i 0.8078 + 0.9877i Poles 0.6119 -0.2226 - 0.4721i -0.2226 + 0.4721i Direct form -3 e) Pole-zero plotting from system function (z-domain) Given that Enter the numerator coefficients of X(z) in matrix form: [1 0 -1/4 0 0 1/5] Enter the denominator coefficients of X(z) in matrix form: [1 0 2 0 0 0 -1/2] P o le -Z e r o p lo t o f X (z )

1

I m a g in a ry P a rt

0 .5

0

-0 .5

-1

-1 .5


-1

-0 .5

0 R e a l P a rt

0 .5

1

1 .5

Page 44

08.508 Expt. No. 5

Digital Signal Processing Lab STUDY OF AM, FM AND PWM WAVEFORM AND RANDOM SEQUENCE GENERATOR USING MATLAB

Objectives: a) To generate AM signal without using built in function in MATLAB b) To generate FM signal without using built in functions in MATLAB c) To generate PWM signals without using built in functions in MATLAB d) To generate random sequences using MATLAB Equipments required: PC installed with MATLAB Theory: Amplitude Modulation (AM) signal can be generated by super-imposing a low frequency signal (modulating signal) into a high frequency signal (carrier signal). Let modulating signal is ) and carrier signal is V , where fm is frequency of modulating signal, is peak m V amplitude of modulating signal, f c is frequency of carrier signal and c is peak amplitude of carrier signal. In general form, AM signal is where m is modulation index of FM and it ranges from 0 to 1. Frequency Modulation (FM) signal can be generated by varying frequency of a high frequency signal (carrier signal) according to the instantaneous amplitude of a low frequency signal (modulating signal). General for of FM signal is where mf is the modulation index. AM and FM signals:


Page 45

08.508


AM spectrum:

FM spectrum:

Pulse Width Modulation (PWM) signal can be generated by varying the duty cycle of carrier rectangular pulse according to the average value of the modulating signal. PWM signals:


Page 46

08.508


PWM spectrum:

Random sequence is a sequence of numbers or symbols that lack any pattern i.e. appear random. Random signal is a noise; hence we cannot predict its value. We can only make statements about the probability of it taking a particular value, or range of values. The Probability density function (pdf) p(x) of a random signal or random variable x is defined to be the probability that the random variable x takes a value between x 0 and x0 + δx. We write this as follows:

Some distributions that describe the probability of a random variable taking certain values are: Uniform distribution is a probability distribution whereby finite number of equally spaced outcomes equally likely to happen. If a random variable has any of n possible values k1, k2,……,kn that are equally spaced and equally probable, then it has a discrete uniform distribution. The probability of any outcome ki is 1 / n. Some properties of uniform distribution are: • The parameters of the uniform distribution are the minimum a and maximum b value. The standard uniform distribution has minimum a = 0 and maximum b = 1 • The probability density function for the uniform distribution is

•

The theoretical mean of the uniform distribution is


Page 47

08.508


The variance of the uniform distribution is Normal (or Gaussian) Distribution is a distribution that is often used as a first approximation to describe real-valued random variables that tend to cluster around a single mean value. The graph of the associated probability density function is "bell"-shaped, and is known as the Gaussian function or bell curve. Some properties of normal distributions are: • The parameters for the normal distribution are the mean µ and the variance σ 2. The standard normal distribution has mean µ = 0 and variance σ2 = 1. • The probability density function for the normal distribution is

• The theoretical mean of the normal distribution is µ • The variance of the normal distribution is σ2 Rayleigh Distribution is often observed when the overall magnitude of a vector is related to its directional components. Some properties of the distribution are:


Page 48

08.508


Algorithms: a) AM wave generation and its spectrum 1. Start 2. Enter modulating signal frequency , carrier signal frequency , peak amplitude of modulating signal , modulation index (or peak amplitude of carrier signal), length of the sequences N and sampling rate Fs 3. Compute index n according to length of sequence 4. Compute , and 5. 6. Display waveforms 7. Obtain frequency spectrum of and index for frequency from 0 to fc+fm 8. Display the spectrum 9. Stop b) FM wave generation and its spectrum 1. Start 2. Enter modulating signal frequency , carrier signal frequency , peak amplitude of modulating signal , modulation index (or peak amplitude of carrier signal), length of the sequences N and sampling rate Fs 3. Compute index n according to length of sequence 4. Compute 5. Display waveforms Department of ECE, VKCET

(choose proportionality constant, kf =1000),

,

and Page 49

08.508

Digital Signal Processing Lab 6. Obtain frequency spectrum of frequency from 0 to fc+fm 7. Display the spectrum

and index for

c) PWM wave generation and its spectrum 1. Start 2. Enter modulating signal frequency , carrier signal frequency , peak amplitude of modulating signal length of the sequences N and sampling rate Fs 3. Compute index n according to length of sequence 4. Compute , 5. Generate PWM signal by comparing instantaneous values of and , if vm(n) > vs(n) , PWM signal vp(n) is high, else low 6. Display waveforms 7. Obtain frequency spectrum of PWM with frequency index 8. Display the spectrum Note: MATLAB function SAWTOOTH: Saw tooth and triangle wave generation. SAWTOOTH(T) generates a saw tooth wave with period 2*pi for the elements of time vector T. SAWTOOTH (T) is like SIN(T), only it creates a saw tooth wave with peaks of +1 to -1 instead of a sine wave. SAWTOOTH (T,WIDTH) generates a modified triangle wave where WIDTH, a scalar parameter between 0 and 1, determines the fraction between 0 and 2*pi at which the maximum occurs. The function increases from -1to 1 on the interval 0 to WIDTH*2*pi, then decreases linearly from 1back to -1 on the interval WIDTH*2*pi to 2*pi. Thus WIDTH = .5 gives you a triangle wave, symmetric about time instant pi with peak amplitude of one. SAWTOOTH (T,1) is equivalent to SAWTOOTH(T). d) Uniform random sequence generator 1. Start 2. Enter the length of the sequence and range (max and min) of the distribution 3. Generate the sequence using rand( ) as seed or a user defined seed 4. Display the sequence 5. Test the sequence using hist ( ), bar( ) 6. Compute mean and median 7. Compute pdf

8. Display pdf along with mean and median 9. Stop Note: 1) MATLAB function RAND uniformly distributed pseudo-random numbers. R = RAND(N) returns an N-by-N matrix containing pseudo-random values drawn from a uniform distribution on the unit interval. RAND(M,N)or RAND([M,N]) returns an M-by-N matrix. RAND(M,N,P,...) or RAND([M,N,P,...]) returns an M-by-N-by-P-by-... array. Department of ECE, VKCET

Page 50

08.508

Digital Signal Processing Lab RAND with no arguments returns a scalar. RAND(SIZE(A)) returns an array the same size as A. 2) MATLAB function hist( ) HIST Histogram. N = HIST(Y) bins the elements of Y into 10 equally spaced containers and returns the number of elements in each container. If Y is a matrix, HIST works down the columns. N = HIST(Y,M), where M is a scalar, uses M bins. N = HIST(Y,X), where X is a vector, returns the distribution of Y among bins with centers specified by X. The first bin includes data between -inf and the first center and the last bin includes data between the last bin and inf. Note: Use HISTC if it is more natural to specify bin edges instead. [N,X] = HIST(...) also returns the position of the bin centers in X. [n,y]=hist(x,N); The function hist divides the interval [min(x), max(x)] into N bins and yields the output [n,y], where n is a vector whose elements are the number of samples in each bin, and y is a vector whose elements are the centers of the bins. When used in this manner, the hist function does not produce a graph; instead we use the bar function. The bar function produces a bar graph which for each value of y, there is a bar whose height is proportional to n. 3) MATLAB function bar( ) BAR Bar graph. BAR(X,Y) draws the columns of the M-by-N matrix Y as M groups of N vertical bars. The vector X must not have duplicate values. e) Normal (Gaussian) random sequence generator 1. Start 2. Enter the length of the sequence 3. Enter mean µ and standard deviation σ 4. Generate the sequence using randn( ) as seed or a user defined seed 5. Display the sequence 6. Test the sequence using hist ( ), bar( ) and histfit( ) 7. Compute mean and median 8. Compute pdf

9. Display pdf along with mean and median 10. Stop Note: 1)MATLAB function randn( ) RANDN Normally distributed random numbers. R = RANDN(N) returns an N-by-N matrix containing pseudo-random values drawn from a normal distribution with mean zero and standard deviation one. RANDN(M,N) or RANDN([M,N]) returns an M-by-N matrix. RANDN(M,N,P,...) or RANDN([M,N,P,...]) returns an M-by-N-by-P-by-... array. RANDN with no arguments returns a scalar. RANDN(SIZE(A)) returns an array the same size as A. Department of ECE, VKCET

Page 51

08.508

Digital Signal Processing Lab f) Rayleigh random sequence generator 1. Start 2. Enter the length of the sequence 3. Enter the scale parameter b 4. Generate the sequence using raylrnd( ) as seed or a user defined seed 5. Display the sequence 6. Test the sequence using hist ( ), bar( ) and histfit( ) 7. Compute mean and median 8. Compute pdf

9. Display pdf along with mean and median 10. Stop Note: MATLAB function raylrnd( ) RAYLRND Random arrays from the Rayleigh distribution. R = RAYLRND(B) returns an array of random numbers chosen from the Rayleigh distribution with parameter B. The size of R is the size of B. Programs: a) AM wave generation and its spectrum %AM wave and spectrum clc; clear all close all; %Input N=input('Enter the length of the waves (Only even number): '); Fs=input('Enter sampling rate in samples/sec: '); m=input('Enter the modulation index in % : '); m=m/100; Vm=input('Enter the peak amplitude of modulating signal in V: '); fm=input('Enter the frequency of modulating signal in Hz: '); fc=input('Enter the frequency of carrier signal in Hz: '); %Signals n=0:1:N-1; vmn=Vm*cos(2*pi*fm*n/Fs); vcn=(Vm/m)*cos(2*pi*fc*n/Fs); %Modulation vamn=(Vm/m).*(1+m.*vmn).*cos(2*pi*fc*n/Fs); %Waveform display subplot(4,1,1); stem(n,vmn); xlabel('n'); ylabel('vm(n)'); title('Modulating signal'); subplot(4,1,2); stem(n,vcn); xlabel('n'); Department of ECE, VKCET

Page 52

08.508


ylabel('vc(n)'); title('Carrier signal'); subplot(4,1,3); stem(n,vamn); xlabel('n'); ylabel('vam(n)'); title('AM (DSB-FC) signal'); %Spectrum display f=(0: (N-1)/2)*Fs/N; Vamf=abs(fftshift(fft(vamn,N)))/N; Vamf=Vamf((N/2+1):N); subplot(4,1,4); stem(f,Vamf); xlabel('f'); ylabel('V'); title('AM spectrum'); b) FM wave generation and its spectrum %FM wave and spectrum clc; clear all close all; %Input N=input('Enter the length of the waves: '); Fs=input('Enter sampling rate in samples per second: '); mf=input('Enter the modulation index: '); fm=input('Enter the frequency of modulating signal in Hz: '); fc=input('Enter the frequency of carrier signal in Hz: '); Vc=input('Enter the peak amplitude of carrier signal in V: '); %Signals n=0:1:N-1; kf=1000; %proprtionality constant of modulator delta=mf*fm; Vm=delta/kf; vmn=Vm*cos(2*pi*fm*n/Fs); vcn=Vc*cos(2*pi*fc*n/Fs); %Modulation vfmn=Vc*cos((2*pi*fc*n/Fs)+mf*sin(2*pi*fm*n/Fs)); %Waveform display subplot(4,1,1); plot(n,vmn); grid on; xlabel('n'); ylabel('vm(n)'); title('Modulating signal'); subplot(4,1,2); plot(n,vcn); Department of ECE, VKCET

Page 53

08.508


grid on; xlabel('n'); ylabel('vc(n)'); title('Carrier signal'); subplot(4,1,3); plot(n,vfmn); xlabel('n'); ylabel('vfm(n)'); title('FM signal'); %Spectrum display f=[0: (N-1)/2]*Fs/N; Vfmf=abs(fftshift(fft(vfmn,N)))/N; Vfmf=Vfmf((N/2+1):N); subplot(4,1,4); stem(f,Vfmf); xlabel('f'); ylabel('V'); title(‘FM spectrum'); c) PWM wave generation and its spectrum %PWM wave and spectrum clc; clear all close all; %Input N=input('Enter the length of the waves: '); Fs=input('Enter sampling rate in samples per second: '); fm=input('Enter the frequency of modulating signal in Hz: '); Vm=input('Enter the peak amplitude of modulating signal (Vm) in V: '); fc=input('Enter the frequency of carrier signal in Hz: '); Vp=input('Enter the peak amplitude of carrier pulse (Vp>Vm) in V: '); %Signals n=0:1:N-1; vmn=Vm*sin(2*pi*fm*n/Fs); vpn=Vp*sawtooth(2*pi*fc*n/Fs); %Modulation vpwmn=zeros(1,N); for i=1:N if(vpn(i)
Page 54

08.508


ylabel('vm(n)'); title('Modulating signal & Sawtooth wave'); subplot(3,1,2); plot(n,vpwmn); axis([0 N 0 10]); xlabel('n'); ylabel('vpwm(n)'); title('PWM signal'); %Spectrum display f=[0: (N-1)/2]*Fs/N; Vpwmf=abs(fftshift(fft(vpwmn,N)))/N; Vpwmf=Vpwmf((N/2+1):N); subplot(3,1,3); stem(f,Vpwmf); xlabel('f'); ylabel('V'); title('PWM spectrum'); d)Uniform random sequence generator %Uniform random sequence generator clc; clear all; close all; %Input N=input('Enter the length of the sequence: '); a=input('Enter the minimum value of range: '); b=input('Enter the maximum value of range: '); %Random sequence x=a+(b-a)*rand(1,N); %Display n=0:N-1; subplot(3,1,1); stem(n,x); xlabel('n'); ylabel('x(n)'); title('Unifrom random sequence'); %Testing disp('Enter any key to test the sequence'); pause; [n,y]=hist(x); subplot(3,1,2); bar(y,n); title('Histogram'); disp('Enter any key to display pdf'); pause; %Calculate pdf for Uniform distribution and display for i=1:length(y) Department of ECE, VKCET

Page 55

08.508


f(i)=1/(b-a); end y=[0:a, y, b:10]; len_f=length(f); len_y=length(y); f=[zeros(1,a+1),f,zeros(1,len_y-len_f-a-1)]; subplot(3,1,3); plot(y,f); xlabel('x'); ylabel('f(x)'); title('pdf'); %mean & variance m=(a+b)/2; var=(b-a)^2/12; text(N/2,1,['Mean = ',num2str(m)]); text(N/2,0.5,['Variance = ',num2str(var)]); e) Normal random sequence generator %Normal random sequence generator clc; clear all; close all; %Input N=input('Enter the length of the sequence: '); sigma=input('Enter standard deviation: '); m=input('Enter mean: '); %Random sequence x=m+sigma*randn(1,N); %Display n=0:N-1; subplot(3,1,1); stem(n,x); xlabel('n'); ylabel('x(n)'); title('Normal random sequence'); %Testing disp('Enter any key to test the sequence'); pause; [n,y]=hist(x); subplot(3,1,2); bar(y,n); %Histogram xlabel('x'); ylabel('bins'); title('Histogram'); disp('Enter any key to display pdf'); pause; Department of ECE, VKCET

Page 56

08.508


%Calculate pdf for Normal distribution p1=-0.5*((y-m)./sigma).^2; p2=sigma.*sqrt(2*pi); f=(1/p2).*exp(p1); subplot(3,1,3); plot(y,f); xlabel('x'); ylabel('f(x)'); title('pdf'); %mean & variance var =sigma^2; text(m,max(f)/2,['Mean = ',num2str(m)]); text(m,max(f)/3,['Variance = ',num2str(var)]); e) Rayleigh random sequence generator %Rayleigh random sequence generator clc; clear all; close all; %Input N=input('Enter the length of the sequence: '); % Number of points b=input('Enter the scale parameter: '); % or standard deviation var=(2-(pi/2))/b^2; % variance sigma=sqrt(var); % generate Rayleigh random sequence x=raylrnd(b,[1,N]); %display n=0:N-1; subplot(3,1,1); stem(n,x); xlabel('n'); ylabel('x(n)'); title('Rayleigh random sequence'); disp('Enter any key to test the sequence'); pause; %histogram [n,y] = hist(x); subplot(3,1,2); bar(y,n); %Histogram xlabel('x'); ylabel('bins'); title('Histogram'); disp('Enter any key to display pdf'); pause; %Calculate pdf for Rayleigh distribution a1=y./b^2; a2=exp(-y.^2./(2*b^2)); Department of ECE, VKCET

Page 57

08.508


f=a1.*a2; y=[0 y]; f=[0 f]; subplot(3,1,3); plot(y,f) ylabel('f(x)') xlabel('x') title('pdf'); %mean & variance m=sqrt(b*pi/2); text(max(y)/2,0.75,['Mean = ',num2str(m)]); text(max(y)/2,0.4,['Variance = ',num2str(var)]); Results: a) AM wave generation and its spectrum Enter the length of the waves (Only even number): 200 Enter sampling rate in samples/sec: 20000 Enter the modulation index in % : 95 Enter the peak amplitude of modulating signal in V: 1 Enter the frequency of modulating signal in Hz: 200 Enter the frequency of carrier signal in Hz: 2000


Page 58

08.508


b) FM wave generation and its spectrum Enter the length of the waves: 200 Enter sampling rate in samples per second: 20000 Enter the modulation index: .5 Enter the frequency of modulating signal in Hz: 200 Enter the frequency of carrier signal in Hz: 2000 Enter the peak amplitude of carrier signal in V: 2


Page 59

08.508


Enter the length of the waves: 200 Enter sampling rate in samples per second: 20000 Enter the modulation index: 4 Enter the frequency of modulating signal in Hz: 200 Enter the frequency of carrier signal in Hz: 2000 Enter the peak amplitude of carrier signal in V: 2

c) PWM wave generation and its spectrum Enter the length of the waves: 200 Enter sampling rate in samples per second: 20000 Enter the frequency of modulating signal in Hz: 200 Enter the peak amplitude of modulating signal (Vm) in V: 4 Enter the frequency of carrier signal in Hz: 1000 Enter the peak amplitude of carrier pulse (Vp>Vm) in V: 8 Department of ECE, VKCET

Page 60

08.508


d)Uniform random sequence generator Enter the length of the sequence: 50 Enter the minimum value of range: 0 Enter the maximum value of range: 5 Enter any key to test the sequence Enter any key to display pdf

Enter the length of the sequence: 50 Enter the minimum value of range: 1 Enter the maximum value of range: 10 Enter any key to test the sequence Enter any key to display pdf


Page 61

08.508


e) Normal random sequence generator Enter the length of the sequence: 50 Enter standard deviation: .45 Enter mean: 0 Enter any key to test the sequence Enter any key to display pdf

Enter the length of the sequence: 50 Enter standard deviation: 0.707 Enter mean: -2 Enter any key to test the sequence Enter any key to display pdf Department of ECE, VKCET

Page 62

08.508


e) Rayleigh random sequence generator Enter the length of the sequence: 50 Enter the scale parameter: .5 Enter any key to test the sequence Enter any key to display pdf

Enter the length of the sequence: 50 Enter the scale parameter: 2 Enter any key to test the sequence Enter any key to display pdf Department of ECE, VKCET

Page 63

08.508



Page 64

08.508 Expt. No. 6

Digital Signal Processing Lab IIR FILTER DESIGN USING MATLAB

Objectives: a) To design IIR digital filter using Butterworth analog filter b) To design IIR digital filter using Chebyshev (Type I and II) analog filters Equipments required: PC installed with MATLAB Theory: Filters are frequency selective network and are used for signal separation and restoration. Filters can be classified into: Analog filters- Realized by R, C and active devices, they are cheap, fast and large dynamic range, but poor performance. Digital filters- Realized by digital system like memory, adder and multipliers; they have better performance in low frequency as well as high frequency operation. They are two types: IIR (Infinite Impulse Response) and FIR (Finite Impulse Response) Ideal frequency response of filters:

Practical frequency response of filters:


Page 65

08.508


Ideal frequency response of digital filters:

IIR Filter Basics: IIR digital filters which can be realized and characterized by the recursive equation,

Where h[k] is the impulse response of the filter, b k and ak are the filter coefficients of the filter, and x[n] and y[n] are the input and output to the filter. The transfer function for the IIR filter is,

The transfer function of the IIR filter, H(z) can be represented as

Where K is the gain, z1, z2, …….are the zeros and p1, p2, …….are the poles of the transfer function, H(z). IIR filter design through analog filters: IIR filters are generally designed by first designing its analog counterparts and then transforming it to digital filters. This method is preferred as analog filter design procedures like Butterworth and Chebyshev. Department of ECE, VKCET

Page 66

08.508


In case of IIR filter design, the most common practice is to convert the digital filter specifications into analog low pass prototype filter specifications. Then to transform it into the desired digital filter transfer function H(z) by a) Impulse invariant b) Bilinear transformation Impulse Invariant Transformation Impulse invariant technique is used to obtain the digital filter whose impulse response is as close as possible to analog filter's impulse response. In the impulse invariant method, starting with the analog filter transfer function, H(s), the impulse response, h(t) is obtained using the Laplace transform. h(t) is then suitably sampled to produce h(nT), and the desired transfer function H(z) is then obtained by z-transforming h(nT), where T is the sampling interval. By transformation:

Then mapping formulae for the transformation is given by

Relation between Analog and Digital Frequencies: MATLAB function: IMPINVAR Impulse invariance method for analog to digital filter conversion. [BZ,AZ] = IMPINVAR(B,A,Fs) creates a digital filter with numerator and denominator coefficients BZ and AZ respectively whose impulse response is equal to the impulse response of the analog filter with coefficients B and A sampled at a frequency of Fs Hertz. The B and A coefficients will be scaled by 1/Fs. If you don't specify Fs, it defaults to 1 Hz. % EXAMPLE: Illustrate the relationship between digital and analog % magnitude responses. [b,a] = butter(4,0.3,'s'); [bz,az] = impinvar(b,a,10); [Ha,Wa] = freqs(b,a,512); Department of ECE, VKCET

Page 67

08.508


[Hz,Wz] = freqz(bz,az,512,10); plot(Wa/(2*pi),20*log10(abs(Ha))); hold on; plot(Wz,20*log10(abs(Hz))) Bilinear Transformation Bilinear transformation is much more important and useful than Impulse Invariant Transformation method, especially in high pass, band pass and stop filters. The method is derived by approximating the first order differential equation with a difference equation. This transformation maps the analog transfer function, H(s) from the s-plane onto the discrete transfer function, H(z) in the z-plane. The entire jΩ axis in the s-plane is mapped onto the z-plane. The left half s-plane is mapped inside the unit circle and the right half s-plane is mapped outside the unit circle in the z-plane. The mapping formulae is

Relation between Analog and Digital Frequencies: MATLAB functions: BILINEAR Bilinear transformation with optional frequency prewarping. [Zd,Pd,Kd] = BILINEAR(Z,P,K,Fs) converts the s-domain transfer function specified by Z, P, and K to a z-transform discrete equivalent obtained from the bilinear transformation: where column vectors Z and P specify the zeros and poles, scalar K specifies the gain, and Fs is the sample frequency in Hz. [NUMd,DENd] = BILINEAR(NUM,DEN,Fs), where NUM and DEN are row vectors containing numerator and denominator transfer function coefficients, NUM(s)/DEN(s), in descending powers of s, transforms to z-transform coefficients NUMd(z)/DENd(z). Butterworth Approximation The magnitude response of a practical low pass filter

Where δ is gain and ω is frequency in radians/sec Department of ECE, VKCET

Page 68

08.508


The magnitude response of Butterworth low pass filter

The low pass Butterworth filter is characterized by the following magnitude squared magnitude response,

Where N is the order of the filter, Ω is analog frequency in radians; ΩC is 3dB cut-off frequency in radians.

The advantage of this approximation is flat response over stop band and pass band. The disadvantage is wide transition band. The transfer function of the filter is: For N is even

For N is odd


Page 69

08.508

Digital Signal Processing Lab The coefficients bk and ck are given by and ck=1 The parameter Bk can be obtained from: for N is even for N is odd

MATLAB functions: 1) 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. 2) 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,'low') designs a lowpass filter. [B,A] = BUTTER(N,Wn,'stop') is a bandstop filter if Wn = [W1 W2]. When used with three left-hand arguments, as in [Z,P,K] = BUTTER(...), the zeros and poles are returned in length N column vectors Z and P, and the gain in scalar K. 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. Chebyshev Approximation The Butterworth approximation, which is monotonic in both the pass band and stop band, gives us a magnitude response which exceeds the specifications for the pass band Department of ECE, VKCET

Page 70

08.508


ripple. The problem is that the order of the filter will be high. Usually, we can tolerate some amount of ripple, and it would be advantageous to have a trade-off between the ripple and the filter order, so that we can obtain a filter which meets the specifications and has a lower order than its Butterworth counterpart. For this we use the Chebychev approximation. The characteristic of Chebyshev approximation is that it minimizes the error between the idealized and the actual filter over the range of filter allowing some ripples in pass band. • Type I, with equal ripple in the pass band, monotonic in the stop band. • Type II, with equal ripple in the stop band, monotonic in the pass band. Magnitude response of Type-I filter:

Magnitude response of Type-II filters:

Type I Chebyshev Filter is characterized by the magnitude squared response is given by

Where

is a Chebyshev polynomial which exhibits equal ripple in the pass band,

N is the order of the polynomial as well as that of the filter, and ɛ determines the pass Department of ECE, VKCET

Page 71

08.508


band ripple, which in decibels is given by, Chebyshev polynomial is given by,

The expression for filter parameter ɛ is given by

The order of the filter

Assume The transfer function H(s) is similar to Butterworth filter, but the coefficients bk and ck is given by

The parameter yN is given by

The parameter Bk can be obtained from: For N is even

For N is odd


Page 72

08.508


Type II Chebyshev Filter is characterized by the magnitude squared response is given by

Chebyshev polynomial is give as Type-II The filter parameter ɛ is given by

The order of the filter

All other parameters are similar to Type-I filter MATLAB functions: 1) CHEB1ORD Chebyshev Type I filter order selection. [N, Wp] = CHEB1ORD(Wp, Ws, Rp, Rs) returns the order N of the lowest order digital Chebyshev Type I filter that loses no more than Rp dB in the passband and has at least Rs dB of attenuation in the stopband. 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] CHEB1ORD also returns Wp, the Chebyshev natural frequency to use with CHEBY1 to achieve the specifications. [N, Wp] = CHEB1ORD(Wp, Ws, Rp, Rs, 's') does the computation for an analog filter, in which case Wp and Ws are in radians/second. 2) CHEB2ORD Chebyshev Type II filter order selection. [N, Ws] = CHEB2ORD(Wp, Ws, Rp, Rs) returns the order N of the lowest order digital Chebyshev Type II 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] CHEB2ORD also returns Wst, the Chebyshev natural frequency to use with CHEBY2 to achieve the specifications. [N, Ws] = CHEB2ORD(Wp, Ws, Rp, Rs, 's') does the computation for an analog filter, in which case Wp and Ws are in radians/second. Department of ECE, VKCET

Page 73

08.508


3) CHEBY1 Chebyshev Type I digital and analog filter design. [B,A] = CHEBY1(N,R,Wp) designs an Nth order lowpass digital Chebyshev filter with R decibels of peak-to-peak ripple in the passband. CHEBY1 returns the filter coefficients in length N+1 vectors B (numerator) and A (denominator). The passband-edge frequency Wp must be 0.0 < Wp < 1.0, with 1.0 corresponding to half the sample rate. Use R=0.5 as a starting point, if you are unsure about choosing R. If Wp is a two-element vector, Wp = [W1 W2], CHEBY1 returns an order 2N bandpass filter with passband W1 < W < W2. [B,A] = CHEBY1(N,R,Wp,'high') designs a highpass filter. [B,A] = CHEBY1(N,R,Wp,'low') designs a lowpass filter. [B,A] = CHEBY1(N,R,Wp,'stop') is a bandstop filter if Wp = [W1 W2]. When used with three left-hand arguments, as in [Z,P,K] = CHEBY1(...), the zeros and poles are returned in length N column vectors Z and P, and the gain in scalar K. CHEBY1(N,R,Wp,'s'), CHEBY1(N,R,Wp,'high','s') and CHEBY1(N,R,Wp,'stop','s') design analog Chebyshev Type I filters. In this case, Wp is in [rad/s] and it can be greater than 1.0. 4) CHEBY2 Chebyshev Type II digital and analog filter design. [B,A] = CHEBY2(N,R,Wst) designs an Nth order lowpass digital Chebyshev filter with the stopband ripple R decibels down and stopband-edge frequency Wst. CHEBY2 returns the filter coefficients in length N+1vectors B (numerator) and A (denominator). The stopband-edge frequency Wst must be 0.0 < Wst < 1.0, with 1.0 corresponding to half the sample rate. Use R = 20 as a starting point, if you are unsure about choosing R. If Wst is a two-element vector, Wst = [W1 W2], CHEBY2 returns an order 2N bandpass filter with passband W1 < W < W2. [B,A] = CHEBY2(N,R,Wst,'high') designs a highpass filter. [B,A] = CHEBY2(N,R,Wst,'low') designs a lowpass filter. [B,A] = CHEBY2(N,R,Wst,'stop') is a bandstop filter if Wst = [W1 W2]. When used with three left-hand arguments, as in [Z,P,K] = CHEBY2(...), the zeros and poles are returned in length N column vectors Z and P, and the gain in scalar K. CHEBY2(N,R,Wst,'s'), CHEBY2(N,R,Wst,'high','s') and CHEBY2(N,R,Wst,'stop','s') design analog Chebyshev Type II filters. In this case, Wst is in [rad/s] and it can be greater than 1.0. Algorithm: (Note: The following algorithms are common for both approximation and all types of filters) Analog Butterworth approximation 1. Start 2. Enter pass band ripple a1 and stop band ripple a2 in dB 3. Enter pass band edge and stop band edge frequencies (fp and fs) in Hz (or in normalized form in radians) 4. Enter sampling frequency in samples/second 5. Check the entered values are correct 6. Normalize the frequencies and compute analog frequencies (Ωp and Ωs) Department of ECE, VKCET

Page 74

08.508

Digital Signal Processing Lab 7. Find orders and cut-off frequencies of the analog filters using buttord( ) function and display the values for verification 8. Compute the transfer functions in s-domain using butter( ) 9. Transform the transfer functions to z-domain by impinvar( ) and bilinear( ) 10. Compute the frequency responses of each type using freqz( ) function 11. Plot the each magnitude and phase responses 12. Stop

Digital Butterworth approximation 1. Start 2. Enter pass band ripple a1 and stop band ripple a2 in dB 3. Enter pass band edge and stop band edge frequencies (fp and fs) in Hz (or in normalized form in radians) 4. Enter sampling frequency in samples/second 5. Check the entered values are correct 6. Compute normalization to edge frequencies, if the entered frequencies are in Hz. Otherwise go to next step. 7. Find orders and cut-off frequencies of the digital filters using buttord( ) function and display the values for verification 8. Compute the transfer functions in z-domain using butter( ) 9. Compute the frequency responses of each type using freqz( ) function 10. Plot the each magnitude and phase responses 11. Stop Analog Chebyshev approximation 1. Start 2. Enter pass band ripple a1 and stop band ripple a2 in dB 3. Enter pass band edge and stop band edge frequencies (f p and fs) in Hz (or in normalized form in radians) 4. Enter sampling frequency in samples/second 5. Check the entered values are correct 6. Normalize the frequencies and compute analog frequencies (Ωp and Ωs) 7. Find orders and cut-off frequencies of the analog filters using cheb1ord( ) and cheb2ord( ) function and display the values for verification 8. Compute the transfer functions in s-domain using cheby1( ) and cheby2( ) 9. Transform the transfer functions to z-domain by impinvar( ) and bilinear( ) 10. Compute the frequency responses of each type using freqz( ) function 11. Plot the each magnitude and phase responses 12. Stop Digital Chebyshev approximation 1. Start 2. Enter pass band ripple a1 and stop band ripple a2 in dB 3. Enter pass band edge and stop band edge frequencies (f p and fs) in Hz (or in normalized form in radians) Department of ECE, VKCET

Page 75

08.508

Digital Signal Processing Lab 4. Enter sampling frequency in samples/second 5. Check the entered values are correct 6. Compute normalization to edge frequencies, if the entered frequencies are in Hz. Otherwise go to next step. 7. Find orders and cut-off frequencies of the digital filters using cheb1ord( ) and cheb2ord( ) functions and display the values for verification 8. Compute the transfer functions in z-domain using cheby1( ) and cheby2( ) 9. Compute the frequency responses of each type using freqz( ) function 10. Plot the each magnitude and phase responses 11. Stop

Programs: 1.

IIR filter using analog Butterworth filter %Butterworth filter

clc; clear all; close all; disp('IIR filter design using Butterworth Approximation'); %Inputs ftype=input('Enter filter type (LPF, HPF, BPF or BSF): '); switch ftype case 'LPF' fp=input('Enter pass band frequency fp in Hz '); fs=input('Enter stop band frequency fs in Hz '); %Error message if fp<=0 error('Pass band edge must be larger than zero'); end if fs<=fp error('Stop band edge must be larger than pass band edge'); end case 'HPF' fp=input('Enter pass band frequency fp in Hz '); fs=input('Enter stop band frequency fs in Hz '); %Error message if fs<=0 error('Pass band edge must be larger than 0'); end if fs>=fp error('Stop band edge must be smaller than pass band edge'); end case 'BPF' fs=input('Enter stop band frequencies fs1 and fs2 (fs2 >fs1) in Hz '); fp=input('Enter pass band frequencies fp1 and fs2 (fp1 >fs1 and fp2
Page 76

08.508


error('Pass band edge must be larger than stop band edge'); end if fp(2)>=fs(2) error('Pass band edge must be smaller than stop band edge'); end case 'BSF' fs=input('Enter stop band frequencies fs1 and fs2 (fs2 >fs1) in Hz '); fp=input('Enter pass band frequencies fp1 and fs2 (fp1 < fs1 and fp2 > fs2) in Hz '); if (fs(1) || fs(2))<=0 error('Stop band edge must be larger than zero'); end if fs(1)<=fp(1) error('Stop band edge must be larger than pass band edge'); end if fs(2)>=fp(2) error('Stop band edge must be smaller than pass band edge'); end end Fs=input('Enter sampling frequency F in samples/sec '); a1=input('Enter pass band ripple in dB '); a2=input('Enter stop band ripple in dB '); %conversion and normalization of frequencies %pi radians/second wp=2*fp/Fs; ws=2*fs/Fs; %conversion of frequency for transformation trans=input('Enter transformation type (IIT or BT): '); switch trans case 'IIT' %Analog frequencies by IIT Wp=wp*Fs*pi; Ws=ws*Fs*pi; case 'BT' %Analog frequencies by BT Wp=2*Fs*tan(wp*pi/2); Ws=2*Fs*tan(ws*pi/2); end %Order and cut-off frequencies of filter [N,Wc]=buttord(Wp,Ws,a1,a2,'s'); %Analog Filter response switch ftype case 'LPF' [b,a]=butter(N,Wc,'s'); case 'HPF' [b,a]=butter(N,Wc,'high','s'); case 'BPF' [b,a]=butter(N,Wc,'bandpass','s'); Department of ECE, VKCET

Page 77

08.508


case 'BSF' [b,a]=butter(N,Wc,'stop','s'); end %Transformation, Digital filter response & Frequency response switch trans case 'IIT' [bz,az]=impinvar(b,a,Fs); [H,f]=freqz(bz,az,512,Fs); case 'BT' [bz,az]=bilinear(b,a,Fs); [H,f]=freqz(bz,az,512,Fs); end %display N and H(z) disp('Order of the filter is:'); disp(N); disp('Filter coefficients are:'); disp('a'); disp(az); disp('b'); disp(bz); %Plotting responses subplot(2,1,1); plot(f,abs(H)); xlabel('f in Hz'); ylabel('|H(w)|'); title(['Magnitude response of IIR ',ftype,' filter with Butterworth approximation using ',trans]); grid on; subplot(2,1,2); plot(f,angle(H)); xlabel('f in Hz'); ylabel('Phase of H(w)'); title(['Phase response of IIR ',ftype,' filter with Butterworth approximation using ',trans]); grid on; 2.

IIR filter using analog Chebyshev filter %Chebychev filter

clc; clear all; close all; disp('IIR filter design using Chebyshev Approximation'); %Inputs ftype=input('Enter filter type (LPF, HPF, BPF or BSF): '); switch ftype case 'LPF' fp=input('Enter pass band frequency fp in Hz '); fs=input('Enter stop band frequency fs in Hz '); %Error message Department of ECE, VKCET

Page 78

08.508


if fp<=0 error('Pass band edge must be larger than 0'); end if fs<=fp error('Stop band edge must be larger than pass band edge'); end case 'HPF' fp=input('Enter pass band frequency fp in Hz '); fs=input('Enter stop band frequency fs in Hz '); %Error message if fs<=0 error('Pass band edge must be larger than 0'); end if fs>=fp error('Stop band edge must be smaller than pass band edge'); end case 'BPF' fs=input('Enter stop band frequencies fs1 and fs2 (fs2 >fs1) in Hz '); fp=input('Enter pass band frequencies fp1 and fs2 (fp1 >fs1 and fp2=fs(2) error('Pass band edge must be smaller than stop band edge'); end case 'BSF' fs=input('Enter stop band frequencies fs1 and fs2 (fs2 >fs1) in Hz '); fp=input('Enter pass band frequencies fp1 and fs2 (fp1 < fs1 and fp2 > fs2) in Hz '); if (fs(1) || fs(2))<=0 error('Stop band edge must be larger than 0'); end if fs(1)<=fp(1) error('Stop band edge must be larger than pass band edge'); end if fs(2)>=fp(2) error('Stop band edge must be smaller than pass band edge'); end end Fs=input('Enter sampling frequency Fs in samples/sec '); a1=input('Enter pass band ripple in dB '); a2=input('Enter stop band ripple in dB '); %conversion and normalization of frequencies %pi radians/second wp=2*fp/Fs; Department of ECE, VKCET

Page 79

08.508


ws=2*fs/Fs; %conversion of frequency for transformation trans=input('Enter transformation type (IIT or BT): '); switch trans case 'IIT' %Analog frequencies by IIT Wp=wp*Fs*pi; Ws=ws*Fs*pi; case 'BT' %Analog frequencies by BT Wp=2*Fs*tan(wp*pi/2); Ws=2*Fs*tan(ws*pi/2); end %Cheby1 or Cheby2 type=input('Enter Chebyshev filter type (type-1 or type-2): '); switch type case 'type-1' %Order and cut-off frequencies of filter [N,Wc]=cheb1ord(Wp,Ws,a1,a2,'s'); switch ftype case 'LPF' %Analog Filter response [b,a]=cheby1(N,a1,Wc,'s'); case 'HPF' [b,a]=cheby1(N,a1,Wc,'high','s'); case 'BPF' [b,a]=cheby1(N,a1,Wc,'bandpass','s'); case 'BSF' [b,a]=cheby1(N,a1,Wc,'stop','s'); end case 'type-2' %Order and cut-off frequencies of filter [N,Wc]=cheb2ord(Wp,Ws,a1,a2,'s'); switch ftype case 'LPF' %Analog Filter response [b,a]=cheby2(N,a2,Wc,'s'); case 'HPF' [b,a]=cheby2(N,a2,Wc,'high','s'); case 'BPF' [b,a]=cheby2(N,a2,Wc,'bandpass','s'); case 'BSF' [b,a]=cheby2(N,a2,Wc,'stop','s'); end end %Transformation, Digital filter response & Frequency response switch trans case 'IIT' Department of ECE, VKCET

Page 80

08.508


[bz,az]=impinvar(b,a,Fs); [H,f]=freqz(bz,az,512,Fs); case 'BT' [bz,az]=bilinear(b,a,Fs); [H,f]=freqz(bz,az,512,Fs); end %display N and H(z) disp('Order of the filter is:'); disp(N); disp('Filter coefficients are:'); disp('a'); disp(az); disp('b'); disp(bz); %Plotting responses subplot(2,1,1); plot(f,20*log10(abs(H))); % Gain in dB xlabel('f in Hz'); ylabel('|H(w)|'); title(['Magnitude response of IIR ',ftype,' filter with Chebyshev ',type,' approximation using ',trans]); grid on; subplot(2,1,2); plot(f,angle(H)); xlabel('f in Hz'); ylabel('Phase of H(w)'); title(['Phase response of IIR ',ftype,' filter with Chebyshev ',type,' approximation using ',trans]); Results: 1. IIR low pass filter using analog Butterworth filter with Impulse Invariant transformation IIR filter design using Butterworth Approximation Enter filter type (LPF, HPF, BPF or BSF): 'LPF' Enter pass band frequency fp in Hz 200 Enter stop band frequency fs in Hz 800 Enter sampling frequency F in samples/sec 4000 Enter pass band ripple in dB . 5 Enter stop band ripple in dB 30 Enter transformation type (IIT or BT): 'IIT' Order of the filter is: 4 Filter coefficients are: a 1.0000 -2.6416 2.7721 -1.3415 0.2503 b 0.0000 0.0092 0.0257 0.0046 0


Page 81

08.508


2. IIR low pass filter using analog Butterworth filter with Bilinear transformation IR filter design using Butterworth Approximation Enter filter type (LPF, HPF, BPF or BSF): 'LPF' Enter pass band frequency fp in Hz 200 Enter stop band frequency fs in Hz 800 Enter sampling frequency F in samples/sec 4000 Enter pass band ripple in dB .5 Enter stop band ripple in dB 30 Enter transformation type (IIT or BT): 'BT' Order of the filter is: 3 Filter coefficients are: a 1.0000 -2.1033 1.5667 -0.4019 b 0.0077 0.0231 0.0231 0.0077

3. IIR highpass filter using analog Butterworth filter with bilinear transformation IIR filter design using Butterworth Approximation Enter filter type (LPF, HPF, BPF or BSF): 'HPF' Enter pass band frequency fp in Hz 1800 Enter stop band frequency fs in Hz 800 Enter sampling frequency F in samples/sec 5000 Enter pass band ripple in dB .25 Enter stop band ripple in dB 20 Enter transformation type (IIT or BT): 'BT' Order of the filter is: 3 Filter coefficients are: a 1.0000 0.3060 0.3581 0.0284 b 0.1280 -0.3839 0.3839 -0.1280


Page 82

08.508


4. IIR band pass filter using analog Butterworth filter using Bilinear transformation IIR filter design using Butterworth Approximation Enter filter type (LPF, HPF, BPF or BSF): 'BPF' Enter stop band frequencies fs1 and fs2 (fs2 >fs1) in Hz [300 5000] Enter pass band frequencies fp1 and fs2 (fp1 >fs1 and fp2fs1) in Hz [4000 6000] Enter pass band frequencies fp1 and fs2 (fp1 < fs1 and fp2 > fs2) in Hz [2000 8000] Enter sampling frequency F in samples/sec 20000 Enter pass band ripple in dB 1 Enter stop band ripple in dB 60 Enter transformation type (IIT or BT): 'BT' Order of the filter is: 6 Filter coefficients are: a 1.0000 -0.0000 -0.1024 -0.0000 0.7816 -0.0000 -0.0484 -0.0000 0.1152 -0.0000 -0.0032 -0.0000 0.0018 b 0.0273 -0.0000 0.1635 -0.0000 0.4089 -0.0000 0.5452 -0.0000 0.4089 -0.0000 0.1635 -0.0000 0.0273 Department of ECE, VKCET

Page 83

08.508


6. IIR low pass filter using analog Chebyshev Type 1 filter with Impulse Invariant transformation IIR filter design using Chebyshev Approximation Enter filter type (LPF, HPF, BPF or BSF): 'LPF' Enter pass band frequency fp in Hz 200 Enter stop band frequency fs in Hz 400 Enter sampling frequency F in samples/sec 4000 Enter pass band ripple in dB 10 Enter stop band ripple in dB 30 Enter transformation type (IIT or BT): 'IIT' Enter Chebyshev filter type (type-1 or type-2): 'type-1' Order of the filter is: 3 Filter coefficients are: a 1.0000 -2.8590 2.7951 -0.9336 b 0 0.0013 0.0012 0 7. IIR low pass filter using analog Chebyshev Type 1 filter with Bilinear transformation IIR filter design using Chebyshev Approximation Enter filter type (LPF, HPF, BPF or BSF): 'LPF' Enter pass band frequency fp in Hz 200 Enter stop band frequency fs in Hz 400 Enter sampling frequency F in samples/sec 4000 Enter pass band ripple in dB 10 Enter stop band ripple in dB 30 Enter transformation type (IIT or BT): 'BT' Enter Chebyshev filter type (type-1 or type-2): 'type-1' Order of the filter is: 3 Filter coefficients are: a 1.0000 -2.8588 2.7949 -0.9337 b 1.0e-003 * [0.3140 0.9421 0.9421 0.31408]


Page 84

08.508


6. IIR highpass filter using analog Chebyshev Type 1 filter with bilinear transformation IIR filter design using Chebyshev Approximation Enter filter type (LPF, HPF, BPF or BSF): 'HPF' Enter pass band frequency fp in Hz 4000 Enter stop band frequency fs in Hz 3500 Enter sampling frequency F in samples/sec 10000 Enter pass band ripple in dB 5 Enter stop band ripple in dB 40 Enter transformation type (IIT or BT): 'BT' Enter Chebyshev filter type (type-1 or type-2): 'type-1' Order of the filter is: 5 Filter coefficients are: a 1.0000 4.2688 7.7135 7.3367 3.6676 0.7718 b 0.0001 -0.0006 0.0012 -0.0012 0.0006 -0.0001 7. IIR band pass filter using analog Chebyshev Type 1 filter using Bilinear transformation IIR filter design using Chebyshev Approximation Enter filter type (LPF, HPF, BPF or BSF): 'BPF' Enter stop band frequencies fs1 and fs2 (fs2 >fs1) in Hz [4000 7000] Enter pass band frequencies fp1 and fs2 (fp1 >fs1 and fp2

Page 85

08.508


8. IIR band stop filter using analog Chebyshev Type 1 filter using Bilinear transformation IIR filter design using Chebyshev Approximation Enter filter type (LPF, HPF, BPF or BSF): 'BSF' Enter stop band frequencies fs1 and fs2 (fs2 >fs1) in Hz [4000 6000] Enter pass band frequencies fp1 and fs2 (fp1 < fs1 and fp2 > fs2) in Hz [2000 8000] Enter sampling frequency F in samples/sec 25000 Enter pass band ripple in dB 10 Enter stop band ripple in dB 50 Enter transformation type (IIT or BT): 'BT' Enter Chebyshev filter type (type-1 or type-2): 'type-1' Order of the filter is: 3 Filter coefficients are: a 1.0000 -0.7481 -0.4019 -0.1062 0.6874 0.3918 -0.6753 b 0.0584 -0.1118 0.2466 -0.2388 0.2466 -0.1118 0.0584 9. IIR low pass filter using analog Chebyshev Type 2 filter using Impulse Invariant Transformation IIR filter design using Chebyshev Approximation Enter filter type (LPF, HPF, BPF or BSF): 'LPF' Enter pass band frequency fp in Hz 3000 Enter stop band frequency fs in Hz 4000 Enter sampling frequency F in samples/sec 12000 Enter pass band ripple in dB 1 Enter stop band ripple in dB 40 Enter transformation type (IIT or BT): 'IIT' Enter Chebyshev filter type (type-1 or type-2): 'type-2' Order of the filter is: 8 Filter coefficients are: a 1.0000 0.2791 1.0094 0.1278 0.2229 -0.0027 0.0090 -0.0008 0.0001 b -0.0987 0.1967 0.4272 0.7601 0.6158 0.3262 0.0721 0.0039 0.0000 Department of ECE, VKCET

Page 86

08.508


10. IIR low pass filter using analog Chebyshev Type 2 filter using Bilinear Transformation IIR filter design using Chebyshev Approximation Enter filter type (LPF, HPF, BPF or BSF): 'LPF' Enter pass band frequency fp in Hz 3500 Enter stop band frequency fs in Hz 4000 Enter sampling frequency F in samples/sec 12000 Enter pass band ripple in dB 10 Enter stop band ripple in dB 30 Enter transformation type (IIT or BT): 'BT' Enter Chebyshev filter type (type-1 or type-2): 'type-2' Order of the filter is: 4 Filter coefficients are: a 1.0000 0.4405 0.6576 0.1669 0.0432 b 0.1964 0.5684 0.7787 0.5684 0.1964 11. IIR highpass filter using analog Chebyshev Type 2 filter with bilinear transformation IIR filter design using Chebyshev Approximation Enter filter type (LPF, HPF, BPF or BSF): 'HPF' Enter pass band frequency fp in Hz 3000 Enter stop band frequency fs in Hz 2500 Enter sampling frequency F in samples/sec 10000 Enter pass band ripple in dB 1 Enter stop band ripple in dB 60 Enter transformation type (IIT or BT): 'BT' Enter Chebyshev filter type (type-1 or type-2): 'type-2' Order of the filter is: 10 Filter coefficients are: a 1.0000 0.7246 1.8775 0.6880 1.0411 0.1277 0.2197 -0.0112 0.0180 -0.0021 0.0003 b 0.0174 -0.0707 0.1806 -0.3280 0.4609 -0.5145 0.4609 -0.3280 0.1806 -0.0707 0.0174 Department of ECE, VKCET

Page 87

08.508


12. IIR band pass filter using analog Chebyshev Type 2 filter using Bilinear transformation IIR filter design using Chebyshev Approximation Enter filter type (LPF, HPF, BPF or BSF): 'BPF' Enter stop band frequencies fs1 and fs2 (fs2 >fs1) in Hz [3000 7000] Enter pass band frequencies fp1 and fs2 (fp1 >fs1 and fp2fs1) in Hz [5000 6000] Enter pass band frequencies fp1 and fs2 (fp1 < fs1 and fp2 > fs2) in Hz [4000 7000] Enter sampling frequency F in samples/sec 20000 Enter pass band ripple in dB 10 Enter stop band ripple in dB 50 Enter transformation type (IIT or BT): 'BT' Enter Chebyshev filter type (type-1 or type-2): 'type-2' Order of the filter is: 3 Filter coefficients are: a 1.0000 0.5813 0.7904 0.3674 0.5075 0.1068 0.0706 b 0.2807 0.2634 0.9036 0.5287 0.9036 0.2634 0.2807 Department of ECE, VKCET

Page 88

08.508 Expt. No. 7

Digital Signal Processing Lab FIR FILTER DESIGN USING MATLAB

Objectives: a) To design FIR digital filter using window methods b) To design FIR digital filter using frequency sampling method Equipments required: PC installed with MATLAB Theory: FIR Filter Basics: The difference equation of FIR filter is

Advantages of FIR filter in design and implementation: • The phase response can be exactly linear. • They are relatively easy to design since there are no stability problem • They are efficient to implement. • The DFT can be used in their implementation. Advantages of a linear-phase response are: • Design problem contains only real arithmetic and not complex arithmetic. • Linear-phase filters provide no delay distortion and only a 6x4 amount of delay • For the filter of length M (or order N=M - 1) the number of operations are of the order of M/2 as in the linear-phase filter implementation. 1. Properties of Linear-Phase FIR filters (Symmetric and Anti-symmetric impulse response) The impulse response of filter h(n), 0 ≤ n ≤ M-1 of length M. The system function of the filter is

Which has M-1 poles at origin z=0 and M-1 zeros at anywhere in the z-plane. The frequency response of the filter is The condition for

linear-phase

Where α is a constant phase delay and h(n) is symmetric about α. Department of ECE, VKCET

Page 89

08.508


Another type,

The impulse response h(n) is anti-symmetric. 2. Linear-phase FIR filter design using Windowing techniques The basic idea behind the window design is to choose a proper ideal frequency-selective filter (which always has a non-causal, infinite-duration impulse response) and then truncate (or window) its impulse response to obtain a linear-phase and causal FIR filter. Ideal frequency selective filter by has unity magnitude gain and linear-phase characteristics over its pass band, and zero response over its stop band. An ideal LPF of bandwidth ωc < π is given by

Where ωc is cut-off frequency, and α is called the sample delay. The impulse response of this filter is of infinite duration and is given by

hd(n) is symmetric with respect to α but not causal. For that truncation (windowing) operation is used. To obtain a causal and linear-phase FIR filter h(n) of length M, we must have

In

general

where w(n) is called window function

In the frequency domain the causal FIR filter response of convolution of and the window response W(ejω) , that is


is given by the periodic

Page 90

08.508


The pictorial representation of the windowing in frequency domain is shown below.

Observations for the above operation are: 1. Since the window w(n) has a finite length equal to M, its response has a peaky main lobe whose width is proportional to

and has side lobes of smaller heights.

2. The periodic convolution produces a smeared version of the ideal response . 3. The main lobe produces a transition band in whose width is responsible for the transition width. This width is then proportional to

. The wider the main lobe, the wider

will be the transition width. 4. The side lobes produce ripples that have similar shapes in both the pass band and stop band. 3. FIR filter design through windows: 1. Rectangular window

There are problems in rectangular windowing method: i) Minimum stop band attenuation is 21dB ii) Due to direct truncation, desired impulse response suffers from Department of ECE, VKCET

Page 91

08.508


Gibb’s phenomenon. That is if we increase M, the width of each side lobe will decrease, but the area under each lobe will remain constant. Therefore the relative amplitudes of side lobes will remain constant, and the minimum stop band attenuation will remain at 21 dB. This implies that all ripples will bunch up near the band edges. Gibb’s phenomenon is shown below.

The other types windows are used to reduce this problems. 2. Bartlett (triangular) window

3. Hanning window


Page 92

08.508


4. Hamming window

5. Blackman window


Page 93

08.508


6. Kaiser window

Where Io[.] is the modified zero order Bessel function and β is a parameter that depends on M and that can be chosen according to the required transition width and optimum stop band attenuation. Design equations for Kaiser Window are:


Page 94

08.508


MATLAB functions: 1) RECTWIN Rectangular window. W = RECTWIN(N) returns the N-point rectangular window. 2) BARTLETT Bartlett window. W = BARTLETT(N) returns the N-point Bartlett window. 3) HANN Hann window. HANN(N) returns the N-point symmetric Hann window in a column vector. 4) HAMMING Hamming window. HAMMING(N) returns the N-point symmetric Hamming window in a column vector. 5) BLACKMAN Blackman window. BLACKMAN(N) returns the N-point symmetric Blackman window in a column vector. 6) KAISER Kaiser window. W = KAISER(N) returns an N-point Kaiser window in the column vector W. W = KAISER(N,BTA) returns the BETA-valued N-point Kaiser window. If omitted, BTA is set to 0.500. Steps to design FIR filter using windows: 1. Calculate order of the filter N using given specifications pass band ripple R p in dB, stop band ripple As in dB, pass band edge frequency in Hz, stop band edge frequency in Hz and sampling rate in samples/sec. Or Specify the order, sampling rate and cut-off frequencies of the filter 2. Choose the window function and compute its coefficients 3. Compute the ideal filter coefficients according to the order of the filter 4. Compute FIR filter coefficients according to the obtained window function and ideal filter coefficients. 5. If the resulting filter has too wide or too narrow transition region, it is necessary to change the filter order by increasing or decreasing it according to needs. 4. FIR filter design using frequency sampling method: The frequency-sampling method for FIR filter design is the simplest and most direct technique imaginable when a desired frequency response has been specified. In this method, the desired frequency response is provided as in the previous method. Then the given frequency response is sampled at a set of equally spaced frequencies to obtain N samples. Sampling the continuous frequency response Hd(ω) at N points gives us the Department of ECE, VKCET

Page 95

08.508


N-point DFT of Hd(k). Thus by using the IDFT formula, the filter coefficients can be calculated using the following formula

Two design techniques: Type-I design Samples are taken at the frequency Desired frequency response at these frequencies are given by The filter coefficients

For M is odd For M is even

Then filter coefficients are

Type-II design Samples are taken at the frequency For M is odd

For M is even


Page 96

08.508


Then filter coefficients are

Steps to design FIR filter using windows: 1. Calculate order of the filter N using given specifications pass band ripple R p in dB, stop band ripple As in dB, pass band edge frequency in Hz, stop band edge frequency in Hz and sampling rate in samples/sec. Or 2. Specify the order, sampling rate and cut-off frequencies of the filter 3. Compute the ideal filter coefficients according to the specified parameter of the filter 4. Compute sampled frequency response by DFT 5. Compute practical filter coefficients by IDFT Algorithm: a) FIR filter using common windows 1. Start 2. Enter cut-off frequencies, sampling rate and order of the filter (or pass band edge frequency, stop band edge frequency, pass band ripple, stop band ripple and sampling rate. The compute cut-off frequency by ω c= (ωp + ωs)/2 , order of the filter by

where c is the constant and value depends on the

width of main lobe of windows) 3. Compute the length of the filter, M 4. Find coefficients of window sequence w(n) by the functions rectwin( ), bartlett( ), hann( ), hamming( ) and blackman( ) 5. Find desired filter coefficients h(n) by the function fir1( ) 6. Find filter response by the function freqz( ) 7. Plot the window, magnitude response and phase response 8. Stop b) FIR filter using Kaiser Window 1. Start 2. Enter cut-off frequency, sampling rate and order of the filter and parameter β 3. Compute the length of the filter, M 4. Find coefficients of window sequence w(n) by the functions kaiser( ) 5. Find desired filter coefficients h(n) by the function fir1( ) 6. Find filter response by the function freqz( ) Department of ECE, VKCET

Page 97

08.508

Digital Signal Processing Lab Plot the window, magnitude response and phase response 8. Stop 7.

Note: Change ‘ftype’ as ‘high’, ‘stop’ or ‘bandpass’ to the function fir1( ) in step 5 for other type filters MATLAB functions: FIR1 FIR filter design using the window method. B = FIR1(N,Wn) designs an N'th order low pass 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 B = FIR1(N,Wn,WIN) designs an N-th order FIR filter using the N+1 length vector WIN to window the impulse response. b = fir1(N,Wn,'ftype',WIN) designs an Nth order FIR filter using N+1 length vector WIN to window the impulse response. where 'ftype' is: 'high' for a high pass filter with cutoff frequency Wn. 'stop' for a band stop filter, if Wn = [w1 w2]. The stop band frequency range is specified by this interval. 'bandpass' for a band pass filter, if Wn = [w1 w2]. The pass band frequency range is specified by this interval. c) FIR filter using frequency sampling 1. Start 2. Enter cut-off frequency, sampling rate and order of the filter 3. Normalize the cut-off frequency in π radians/sec 4. Generate normalized frequency index (W) ranges from 0 to 1, in that cut-off frequencies must repeat two times for narrow transition band for ideal response (for LPF and HPF assign W=[0 wc wc 1]) 5. Generate amplitudes corresponding to ideal response of filter (A) (for LPF assign A=[1 1 0 0] 6. Find desired filter coefficients h(n) by the function fir2( ) 7. Find filter response by the function freqz( ) 8. Plot the ideal and desired magnitude response and phase response 9. Stop MATLAB functions: FIR2 FIR arbitrary shape filter design using the frequency sampling method. B = FIR2(N,F,A) designs an N'th order FIR digital filter with the frequency response specified by vectors F and A, and returns the filter coefficients in length N+1 vector B. Vectors F and A specify the frequency and magnitude breakpoints for the filter such that PLOT(F,A) would show a plot of the desired frequency response.The frequencies in F must be between 0.0 < F < 1.0, with 1.0 corresponding to half the sample rate. They must be in increasing order and start with 0.0 and end with 1.0. The filter B is real, and has linear phase, i.e., symmetric Department of ECE, VKCET

Page 98

08.508

Digital Signal Processing Lab coefficients obeying B(k) = B(N+2-k), k = 1,2,...,N+1. By default FIR2 windows the impulse response with a Hamming window.

Program a) FIR filter using windows %FIR filter using windows clc; clear all; close all; disp('FIR filter design using windows'); %Inputs ftype=input('Enter filter type (LPF, HPF, BPF or BSF): '); switch ftype case 'LPF' fc=input('Enter cut-off frequency fc in Hz '); %Error message if fc<=0 error('Cut-off frequency must be larger than zero'); end case 'HPF' fc=input('Enter cut-off frequency fc in Hz '); %Error message if fc<=0 error('Cut-off frequency must be larger than zero'); end case 'BPF' fc=input('Enter cut-off frequencies fc1 and fc2 (fc2 >fc1) in Hz '); if (fc(1) || fc(2))<=0 error('Cut-off frequencies must be larger than 0'); end if fc(1)<=fc(2) error('Lower cut-off frequency must be smaller than upper cut-off frequency'); end case 'BSF' fc=input('Enter cut-off frequencies fc1 and fc2 (fc2 >fc1) in Hz '); if (fc(1) || fc(2))<=0 error('Cut-off frequencies must be larger than 0'); end if fc(1)<=fc(2) error('Lower cut-off frequency must be smaller than upper cut-off frequency'); end end Fs=input('Enter sampling frequency F in samples/sec '); N=input('Enter the order of the filter '); %length of filter M=N+1; %conversion and normalization of frequencies Department of ECE, VKCET

Page 99

08.508


%pi radians/second wc=2*fc/Fs; %input to window type wintype=input('Enter the window (Rectangular,Bartlett,Hann,Hamming,Blackman, Kaiser): '); switch wintype case 'Rectangular' %Rectangular window w=rectwin(M); case 'Bartlett' %Bartlett window w=bartlett(M); case 'Hann' %Hanning window w=hann(M); case 'Hamming' %Hamming window w=hamming(M); case 'Blackman' %Blackmann window w=blackman(M); case 'Kaiser' %Kaiser window beta=input('Enter shape parameter beta'); w=kaiser(M,beta); end %Digital filter response switch ftype case 'LPF' bz=fir1(N,wc,w); case 'HPF' bz=fir1(N,wc,'high',w); case 'BPF' bz=fir1(N,wc,'bandpass',w); case 'BSF' bz=fir1(N,wc,'stop',w); end %Frequency response [H,f]=freqz(bz,1,512,Fs); %display H(z) disp('Filter coefficients are:'); disp('b'); disp(bz); %Plotting windows and responses n=0:1:M-1; subplot(3,1,1); stem(n,win); Department of ECE, VKCET

Page 100

08.508


xlabel('n'); ylabel('w(n)'); title([wintype,' window for FIR ',ftype]); subplot(3,1,2); plot(f,20*log10(abs(H))); xlabel('f in Hz'); ylabel('H(w)in dB'); title(['Magnitude response of FIR ',ftype,' filter with ',wintype,' window']); grid on; subplot(3,1,3); plot(f,angle(H)); xlabel('f in Hz'); ylabel('Phase of H(w)'); title(['Phase response of FIR ',ftype,' filter with ',wintype,' window']); grid on; b) FIR filter using frequency sampling method %FIR filter using frequency sampling clc; clear all; close all; disp('FIR filter design using frequency sampling'); %Inputs ftype=input('Enter filter type (LPF, HPF, BPF or BSF): '); switch ftype case 'LPF' fc=input('Enter cut-off frequency fc in Hz '); %Error message if fc<=0 error('Cut-off frequency must be larger than zero'); end case 'HPF' fc=input('Enter cut-off frequency fc in Hz '); %Error message if fc<=0 error('Cut-off frequency must be larger than zero'); end case 'BPF' fc=input('Enter cut-off frequencies fc1 and fc2 (fc2 >fc1) in Hz '); if (fc(1) || fc(2))<=0 error('Cut-off frequencies must be larger than 0'); end if fc(1)>=fc(2) error('Lower cut-off frequency must be smaller than upper cut-off frequency'); end case 'BSF' Department of ECE, VKCET

Page 101

08.508

Digital Signal Processing Lab fc=input('Enter cut-off frequencies fc1 and fc2 (fc2 >fc1) in Hz '); if (fc(1) || fc(2))<=0 error('Cut-off frequencies must be larger than 0'); end if fc(1)>=fc(2) error('Lower cut-off frequency must be smaller than upper cut-off frequency'); end

end Fs=input('Enter sampling frequency Fs in samples/sec '); N=input('Enter the order of the filter '); %conversion and normalization of frequencies %pi radians/second wc=2*fc/Fs; switch ftype case 'LPF' % Vector of frequency break points (omega values in units of pi) W=[0 wc wc 1]; % Amplitudes corresponding to Filter A=[1 1 0 0]; case 'HPF' % Vector of frequency break points (omega values in units of pi) W=[0 wc wc 1]; % Amplitudes corresponding to Filter A=[0 0 1 1]; case 'BPF' % Vector of frequency break points (omega values in units of pi) W=[0 wc(1) wc(1) wc(2) wc(2) 1]; % Amplitudes corresponding to Filter A=[0 0 1 1 0 0]; case 'BSF' % Vector of frequency break points (omega values in units of pi) W=[0 wc(1) wc(1) wc(2) wc(2) 1]; % Amplitudes corresponding to Filter A=[1 1 0 0 1 1]; end % Filter coefficients bz=fir2(N,W,A); %Frequency response [H,f]=freqz(bz,1,512,Fs); %display H(z) disp('Filter coefficients are:'); disp('b'); disp(bz); %Plotting desired and practical responses f1=W*Fs/2; subplot(3,1,1); plot(f1,A); Department of ECE, VKCET

Page 102

08.508


axis([0 max(f) 0 2]); xlabel('f'); ylabel('H(w)'); title(['Ideal response of ',ftype,]); grid on; subplot(3,1,2); plot(f,20*log10(abs(H))); xlabel('f in Hz'); ylabel('Hd(w)in dB'); title(['Desired Magnitude response of FIR ',ftype,' using Frequency sampling']); subplot(3,1,3); plot(f,angle(H)); xlabel('f in Hz'); ylabel('Phase of Hd(w)'); title(['Desired Phase response of FIR ',ftype,' using Frequency sampling']); Results: 1. FIR low pass filter using Rectangular window FIR filter design using windows Enter filter type (LPF, HPF, BPF or BSF): ‘LPF’ Enter cut-off frequency fc in Hz 5000 Enter sampling frequency Fs in samples/sec 25000 Enter the order of the filter 20 Enter the window (Rectangular,Bartlett,Hann,Hamming,Blackman,Kaiser): ‘Rectangular’ Filter coefficients are: b -0.0000 -0.0352 -0.0244 0.0279 0.0527 -0.0000 -0.0791 -0.0652 0.0978 0.3164 0.4181 0.3164 0.0978 -0.0652 -0.0791 -0.0000 0.0527 0.0279 -0.0244 -0.0352 -0.0000


Page 103

08.508


2. FIR low pass filter using Bartlett window FIR filter design using windows Enter filter type (LPF, HPF, BPF or BSF): ‘LPF’ Enter cut-off frequency fc in Hz 5000 Enter sampling frequency Fs in samples/sec 25000 Enter the order of the filter 20 Enter the window (Rectangular,Bartlett,Hann,Hamming,Blackman,Kaiser): ‘Bartlett’ Filter coefficients are: b 0 -0.0035 -0.0049 0.0084 0.0211 -0.0000 -0.0475 -0.0456 0.0782 0.2848 0.4181 0.2848 0.0782 -0.0456 -0.0475 -0.0000 0.0211 0.0084 -0.0049 -0.0035 0

3. FIR low pass filter using Hann window FIR filter design using windows Enter filter type (LPF, HPF, BPF or BSF): 'LPF' Enter cut-off frequency fc in Hz 5000 Enter sampling frequency Fs in samples/sec 25000 Enter the order of the filter 20 Enter the window (Rectangular,Bartlett,Hann,Hamming,Blackman,Kaiser): 'Hann' Filter coefficients are: b 0 -0.0008 -0.0022 0.0055 0.0174 -0.0000 -0.0495 -0.0494 0.0845 0.2949 0.3994 0.2949 0.0845 -0.0494 -0.0495 -0.0000 0.0174 0.0055 -0.0022 -0.0008 0 4. FIR low pass filter using Hamming window FIR filter design using windows Enter filter type (LPF, HPF, BPF or BSF): 'LPF' Enter cut-off frequency fc in Hz 5000 Enter sampling frequency Fs in samples/sec 25000 Enter the order of the filter 20 Enter the window (Rectangular,Bartlett,Hann,Hamming,Blackman,Kaiser): 'Hamming' Filter coefficients are: b -0.0000 -0.0035 -0.0039 0.0072 0.0201 -0.0000 -0.0517 -0.0506 0.0855 0.2965 0.4008 0.2965 0.0855 -0.0506 -0.0517 -0.0000 0.0201 0.0072 -0.0039 -0.0035 -0.0000 Department of ECE, VKCET

Page 104

08.508


5. FIR low pass filter using Blackman window FIR filter design using windows Enter filter type (LPF, HPF, BPF or BSF): 'LPF' Enter cut-off frequency fc in Hz 5000 Enter sampling frequency F in samples/sec 25000 Enter the order of the filter 20 Enter the window (Rectangular,Bartlett,Hann,Hamming,Blackman,Kaiser): 'Blackman' Filter coefficients are: b 0.0000 -0.0003 -0.0009 0.0027 0.0101 -0.0000 -0.0386 -0.0430 0.0794 0.2906 0.3999 0.2906 0.0794 -0.0430 -0.0386 -0.0000 0.0101 0.0027 -0.0009 -0.0003 0.0000

6. FIR low pass filter using Kaiser window FIR filter design using windows Enter filter type (LPF, HPF, BPF or BSF): 'LPF' Enter cut-off frequency fc in Hz 5000 Enter sampling frequency F in samples/sec 25000 Enter the order of the filter 20 Enter the window (Rectangular,Bartlett,Hann,Hamming,Blackman,Kaiser): 'Kaiser' Enter shape parameter beta 0 Filter coefficients are: b -0.0000 -0.0352 -0.0244 0.0279 0.0527 -0.0000 -0.0791 -0.0652 0.0978 0.3164 0.4181 0.3164 0.0978 -0.0652 -0.0791 -0.0000 0.0527 0.0279 -0.0244 -0.0352 -0.0000 7. FIR filter design using windows Enter filter type (LPF, HPF, BPF or BSF): 'LPF' Enter cut-off frequency fc in Hz 5000 Enter sampling frequency F in samples/sec 25000 Enter the order of the filter 20 Enter the window (Rectangular,Bartlett,Hann,Hamming,Blackman,Kaiser): 'Kaiser' Enter shape parameter beta 7.5 Filter coefficients are: b -0.0000 -0.0032 -0.0042 0.0077 0.0209 -0.0000 -0.0523 -0.0508 0.0855 0.2962 0.4002 0.2962 0.0855 -0.0508 -0.0523 -0.0000 0.0209 0.0077 -0.0042 -0.0032 -0.0000 Department of ECE, VKCET

Page 105

08.508


8. FIR high pass filter using Kaiser window FIR filter design using windows Enter filter type (LPF, HPF, BPF or BSF): 'HPF' Enter cut-off frequency fc in Hz 5000 Enter sampling frequency F in samples/sec 25000 Enter the order of the filter 20 Enter the window (Rectangular,Bartlett,Hann,Hamming,Blackman,Kaiser): 'Kaiser' Enter shape parameter beta 0 Filter coefficients are: b -0.0000 0.0344 0.0239 -0.0274 -0.0516 0.0000 0.0775 0.0638 -0.0957 -0.3098 0.6141 -0.3098 -0.0957 0.0638 0.0775 0.0000 -0.0516 -0.0274 0.0239 0.0344 -0.0000

9. FIR filter design using windows Enter filter type (LPF, HPF, BPF or BSF): 'HPF' Enter cut-off frequency fc in Hz 5000 Enter sampling frequency F in samples/sec 25000 Enter the order of the filter 20 Enter the window (Rectangular,Bartlett,Hann,Hamming,Blackman,Kaiser): 'Kaiser' Enter shape parameter beta 3 Filter coefficients are: b -0.0000 0.0102 0.0096 -0.0139 -0.0317 0.0000 0.0622 0.0561 -0.0895 -0.3004 0.6027 -0.3004 -0.0895 0.0561 0.0622 0.0000 -0.0317 -0.0139 0.0096 0.0102 -0.0000 10. FIR filter design using windows Enter filter type (LPF, HPF, BPF or BSF): 'HPF' Enter cut-off frequency fc in Hz 5000 Enter sampling frequency F in samples/sec 25000 Enter the order of the filter 20 Enter the window (Rectangular,Bartlett,Hann,Hamming,Blackman,Kaiser): 'Kaiser' Enter shape parameter beta 20 Filter coefficients are: b -0.0000 0.0000 0.0000 -0.0001 -0.0010 0.0000 0.0149 0.0254 -0.0631 -0.2746 0.6000 -0.2746 -0.0631 0.0254 0.0149 0.0000 -0.0010 -0.0001 0.0000 0.0000 -0.0000 Department of ECE, VKCET

Page 106

08.508


11. FIR filter design using windows Enter filter type (LPF, HPF, BPF or BSF): 'BPF' Enter cut-off frequencies fc1 and fc2 (fc2 >fc1) in Hz [4000 8000] Enter sampling frequency Fs in samples/sec 25000 Enter the order of the filter 20 Enter the window (Rectangular,Bartlett,Hann,Hamming,Blackman,Kaiser): 'Kaiser' Enter shape parameter beta 0 Filter coefficients are: b 0.0516 -0.0393 -0.0567 0.0150 -0.0130 0.0244 0.1471 -0.0418 -0.2811 0.0203 0.3374 0.0203 -0.2811 -0.0418 0.1471 0.0244 -0.0130 0.0150 -0.0567 -0.0393 0.0516 12. FIR filter design using windows Enter filter type (LPF, HPF, BPF or BSF): 'BPF' Enter cut-off frequencies fc1 and fc2 (fc2 >fc1) in Hz [4000 8000] Enter sampling frequency F in samples/sec 25000 Enter the order of the filter 20 Enter the window (Rectangular,Bartlett,Hann,Hamming,Blackman,Kaiser): 'Kaiser' Enter shape parameter beta 6 Filter coefficients are: b 0.0007 -0.0020 -0.0066 0.0031 -0.0043 0.0114 0.0906 -0.0316 -0.2445 0.0192 0.3277 0.0192 -0.2445 -0.0316 0.0906 0.0114 -0.0043 0.0031 -0.0066 -0.0020 0.0007 13. FIR filter design using windows Enter filter type (LPF, HPF, BPF or BSF): 'BSF' Enter cut-off frequencies fc1 and fc2 (fc2 >fc1) in Hz [4000 8000] Enter sampling frequency F in samples/sec 25000 Enter the order of the filter 20 Enter the window (Rectangular,Bartlett,Hann,Hamming,Blackman,Kaiser): 'Kaiser' Enter shape parameter beta 1 Filter coefficients are: b -0.0383 0.0305 0.0459 -0.0126 0.0113 -0.0217 -0.1333 0.0385 0.2618 -0.0190 0.6737 -0.0190 0.2618 0.0385 -0.1333 -0.0217 0.0113 -0.0126 0.0459 0.0305 -0.0383 Department of ECE, VKCET

Page 107

08.508


14. FIR filter design using windows Enter filter type (LPF, HPF, BPF or BSF): 'BSF' Enter cut-off frequencies fc1 and fc2 (fc2 >fc1) in Hz [4000 8000] Enter sampling frequency F in samples/sec 25000 Enter the order of the filter 20 Enter the window (Rectangular,Bartlett,Hann,Hamming,Blackman,Kaiser): 'Kaiser' Enter shape parameter beta 9 Filter coefficients are: b -0.0000 0.0004 0.0019 -0.0013 0.0023 -0.0075 -0.0688 0.0269 0.2247 -0.0185 0.6800 -0.0185 0.2247 0.0269 -0.0688 -0.0075 0.0023 -0.0013 0.0019 0.0004 -0.0000

15. FIR filter design using frequency sampling Enter filter type (LPF, HPF, BPF or BSF): 'LPF' Enter cut-off frequency fc in Hz 5000 Enter sampling frequency Fs in samples/sec 25000 Enter the order of the filter 20 Filter coefficients are: b 0.0000 -0.0033 -0.0038 0.0069 0.0197 0.0002 -0.0510 -0.0505 0.0848 0.2958 0.4004 0.2958 0.0848 -0.0505 -0.0510 0.0002 0.0197 0.0069 -0.0038 -0.0033 0.0000

16. FIR filter design using frequency sampling Enter filter type (LPF, HPF, BPF or BSF): 'HPF' Enter cut-off frequency fc in Hz 5000 Enter sampling frequency Fs in samples/sec 25000 Enter the order of the filter 20 Filter coefficients are: b -0.0000 0.0033 0.0038 -0.0069 -0.0197 -0.0002 0.0510 0.0505 -0.0848 -0.2958 0.5996 -0.2958 -0.0848 0.0505 0.0510 -0.0002 -0.0197 -0.0069 0.0038 0.0033 -0.0000


Page 108

08.508


17. FIR filter design using frequency sampling Enter filter type (LPF, HPF, BPF or BSF): 'BPF' Enter cut-off frequencies fc1 and fc2 (fc2 >fc1) in Hz [4000 8000] Enter sampling frequency Fs in samples/sec 25000 Enter the order of the filter 20 Filter coefficients are: b 0.0037 -0.0036 -0.0087 0.0037 -0.0047 0.0120 0.0943 -0.0312 -0.2428 0.0184 0.3203 0.0184 -0.2428 -0.0312 0.0943 0.0120 -0.0047 0.0037 -0.0087 -0.0036 0.0037

18. FIR filter design using frequency sampling Enter filter type (LPF, HPF, BPF or BSF): 'BSF' Enter cut-off frequencies fc1 and fc2 (fc2 >fc1) in Hz [4000 8000] Enter sampling frequency Fs in samples/sec 25000 Enter the order of the filter 20 Filter coefficients are: b -0.0037 0.0036 0.0087 -0.0037 0.0047 -0.0120 -0.0943 0.0312 0.2428 -0.0184 0.6797 -0.0184 0.2428 0.0312 -0.0943 -0.0120 0.0047 -0.0037 0.0087 0.0036 -0.0037


Page 109

08.508 Expt. No. 8

Digital Signal Processing Lab OPTIMUM EQUIRIPPLE FIR FILTER USING MATLAB

Objectives: To design optimum equiripple FIR digital filter using window methods Equipments required: PC installed with MATLAB Theory: The design of FIR filters using Windows methods leads to good performance filters and also easy to understand and implement. However, they have some disadvantages. First, we cannot specify the band frequencies ωp and ωs precisely in the design; that is, we have to accept whatever values we obtain after the design. Second, we cannot specify both δ1 and δ2 ripple factors simultaneously. Either we have δ 1 = δ2 in the window design method, or we can optimize only δ2 in the frequency sampling method. Finally, the approximation error-that is, the difference between the ideal response and the actual response is not uniformly distributed over the band intervals. It is higher near the band edges and smaller in the regions away from band edges. By distributing the error uniformly, we can obtain a lower-order filter satisfying the same specifications. One of the techniques to eliminate these problems is Chebyshev error criterion. The resulting filters are known as Equiripple FIR Filters. FIR filter design algorithm, which uses polynomial interpolation for its solution is known as the Parks McClellan algorithm, and it incorporates the Remez exchange routine for polynomial solution. This algorithm is available as a subroutine on many computing platforms. Mathematical Background Parameters of the four FIR filters types: Type I Order Even F(ω) 1 M N/2 ω0 0

II Odd cos(ω/2) (N-1)/2 0

III Even sin(ω) (N-2)/2 π/2

IV Odd sin(ω/2) (N-1)/2 π/2

In order to minimize the error we need to define an error function E(θ) and a weight function W(θ) which defines the relative importance of the error at any given frequency θ. Then, the error function can be described as follows: where is the desired amplitude response, and is the actual amplitude response. A simple weight function , could be defined as follows:

And the resulting amplitude response

is defined by:

and where F(ω) and M are obtained from the parameters table Department of ECE, VKCET

Page 110

08.508


The problem here is to obtain the coefficients b[k] that minimize the maximum absolute weighted error |E(ω)|, that is, to obtain where is in the operating frequency range of the filter. The Alternation Theorem This theorem states that there exist at least 2 + K frequencies ω i, {0 ≤ i ≤ K+1} where the maximum error, ε occurs. That is, The last equation shows that the sign changes K+1 times, resulting in an oscillation or ripple on the band of interest. The Remez Exchange Algorithm The most common implementation of the Remez Exchange Algorithm is the version by Parks, McClellan and Rabiner. Its objective is to obtain the coefficients b[k] that minimizes ε. It uses the properties of the Alternation Theorem. The first step is to find the order N of the desired filter. The following is empirical formulae proposed by Kaiser:

Where is pass band edge frequency, is stop band edge frequency, deviation, is stop band allowed deviation and

is pass band allowed

Where Rp and As are the attenuations on the pass band and stop band respectively. MATLAB functions: 1. FIRPMORD Parks-McClellan optimal equiripple FIR order estimator. [N,Fo,Ao,W] = FIRPMORD(F,A,DEV,Fs) finds the approximate order N, normalized frequency band edges Fo, frequency band magnitudes Ao and weights W to be used by the FIRPM function as follows: B = FIRPM(N,Fo,Ao,W) The resulting filter will approximately meet the specifications given by the input parameters F, A, and DEV. F is a vector of cutoff frequencies in Hz, in ascending order between 0 and half the sampling frequency Fs. If you do not specify Fs, it defaults to 2. A is a vector specifying the desired function's amplitude on the bands defined by F. The length of F is twice the length of A, minus 2 (it must therefore be even). The first frequency band always starts at zero, and the last always ends at Fs/2. It is not necessary to add these elements to the F vector. DEV is a vector of maximum deviations or ripples (in linear units) allowable for each band. DEV must have the same length as A. EXAMPLE Design a lowpass filter with a passband-edge frequency of 1500Hz, a stopband-edge of Department of ECE, VKCET

Page 111

08.508


2000Hz, passband ripple of 0.01, stopband ripple of 0.1, and a sampling frequency of 8000Hz: [n,fo,mo,w] = firpmord( [1500 2000], [1 0], [0.01 0.1], 8000 ); b = firpm(n,fo,mo,w); This is equivalent to c = firpmord( [1500 2000], [1 0], [0.01 0.1], 8000, 'cell'); b = firpm(c{:}); CAUTION 1: The order N is often underestimated. If the filter does not meet the original specifications, a higher order such as N+1 or N+2 will. CAUTION 2: Results are inaccurate if cutoff frequencies are near zero frequency or the Nyquist frequency. 2. FIRPM Parks-McClellan optimal equiripple FIR filter design. B=FIRPM(N,F,A) returns a length N+1 linear phase (real, symmetric coefficients) FIR filter which has the best approximation to the desired frequency response described by F and A in the minimax sense. F is a vector of frequency band edges in pairs, in ascending order between 0 and 1. 1 corresponds to the Nyquist frequency or half the sampling frequency. At least one frequency band must have a non-zero width. A is a real vector the same size as F which specifies the desired amplitude of the frequency response of the resultant filter B. % Example of a length 31 lowpass filter: h=firpm(30,[0 .1 .2 .5]*2,[1 1 0 0]); % Example of a low-pass differentiator: h=firpm(44,[0 .3 .4 1],[0 .2 0 0],'differentiator'); % Example of a type 4 highpass filter: h=firpm(25,[0 .4 .5 1],[0 0 1 1],'h'); Algorithm: Optimum equiripple FIR filter 1. Start 2. Enter pass band ripple Rp and stop band ripple As in dB 3. Enter pass band edge and stop band edge frequencies (fp and fs) in Hz (or in normalized form in radians) 4. Enter sampling frequency in samples/second 5. Check the entered values are correct 6. Calculate δ1 and δ2 populate vector dev and assign dev=[ δ1 δ2] for LPF, dev=[ δ2 δ1] for HPF, dev=[ δ2 δ1 δ2 ] for BPF and dev=[ δ1 δ2 δ1] for BSF 7. Assign A = [1 0] for LPF, A=[0 1] for HPF, A=[0 1 0] for BPF and A=[1 0 1] for BSF 8. Use the firpmord function to determine order of the filter, normalized frequencies, frequency band magnitudes and weights. 9. Use the firpm function to determine the digital filter coefficients 10. Find filter response by the function freqz( ) 11. Plot the magnitude response and phase response 12. Stop Department of ECE, VKCET

Page 112

08.508


Programs: %FIR filter using Optimum equiripple clc; clear all; close all; disp('FIR filter design using optimum equiripple algorithm'); %Inputs ftype=input('Enter filter type (LPF, HPF, BPF or BSF): '); switch ftype case 'LPF' fp=input('Enter pass band frequency fp in Hz '); fs=input('Enter stop band frequency fs in Hz '); %Error message if fp<=0 error('Pass band edge must be larger than 0'); end if fs<=fp error('Stop band edge must be larger than pass band edge'); end F=[fp fs]; %frequencies in matrix form case 'HPF' fp=input('Enter pass band frequency fp in Hz '); fs=input('Enter stop band frequency fs in Hz '); %Error message if fs<=0 error('Pass band edge must be larger than 0'); end if fs>=fp error('Stop band edge must be smaller than pass band edge'); end F=[fs fp]; %frequencies in matrix form case 'BPF' fs=input('Enter stop band frequencies fs1 and fs2 (fs2 >fs1) in Hz '); fp=input('Enter pass band frequencies fp1 and fs2 (fp1 >fs1 and fp2=fs(2) error('Pass band edge must be smaller than stop band edge'); end F=[fs(1) fp(1) fp(2) fs(2)]; %frequencies in matrix form case 'BSF' fs=input('Enter stop band frequencies fs1 and fs2 (fs2 >fs1) in Hz '); fp=input('Enter pass band frequencies fp1 and fs2 (fp1 < fs1 and fp2 > fs2) in Hz '); Department of ECE, VKCET

Page 113

08.508


if (fs(1) || fs(2))<=0 error('Stop band edge must be larger than 0'); end if fs(1)<=fp(1) error('Stop band edge must be larger than pass band edge'); end if fs(2)>=fp(2) error('Stop band edge must be smaller than pass band edge'); end F=[fp(1) fs(1) fs(2) fp(2)]; %frequencies in matrix form end Fs=input('Enter sampling frequency Fs in samples/sec '); Rp=input('Enter pass band ripple in dB '); As=input('Enter stop band attenuation in dB '); %d1 and d2 calculation d1=(10^(Rp/20)-1)/(10^(Rp/20)+1); d2=10^(-As/20); switch ftype case 'LPF' %populating dev dev=[d1 d2]; % Amplitudes corresponding to Filter A=[1 0]; case 'HPF' %populating dev dev=[d2 d1]; % Amplitudes corresponding to Filter A=[0 1]; case 'BPF' %populating dev dev=[d2 d1 d2]; % Amplitudes corresponding to Filter A=[0 1 0]; case 'BSF' %populating dev dev=[d1 d2 d1]; % Amplitudes corresponding to Filter A=[1 0 1]; end %Filter order, normalized frequency band edges Fo, frequency band %magnitudes Ao and weights W [N,F0,A0,W] = firpmord(F,A,dev,Fs); % Filter coefficients bz=firpm(N,F0,A0,W); %Frequency response [H,f]=freqz(bz,1,512,Fs); %display N and H(z) Department of ECE, VKCET

Page 114

08.508


disp('Filter order: '); disp(N); disp('Filter coefficients are:'); disp('b'); disp(bz); %Plotting desired and practical responses F0=F0*Fs/2; subplot(3,1,1); plot(F0,A0); axis([0 max(f) 0 2]); xlabel('f'); ylabel('H(w)'); title(['Ideal response of ',ftype,]); grid on; subplot(3,1,2); plot(f,20*log10(abs(H))); xlabel('f in Hz'); ylabel('Hd(w)in dB'); title(['Desired Magnitude response of FIR ',ftype,' using Optimum equiripple']); subplot(3,1,3); plot(f,angle(H)); xlabel('f in Hz'); ylabel('Phase of Hd(w)'); title(['Desired Phase response of FIR ',ftype,' using Optimum equiripple']); Result: 1. FIR low pass filter using optimum equiripple FIR filter design using optimum equiripple algorithm Enter filter type (LPF, HPF, BPF or BSF): 'LPF' Enter pass band frequency fp in Hz 4500 Enter stop band frequency fs in Hz 5000 Enter sampling frequency Fs in samples/sec 15000 Enter pass band ripple in dB 1 Enter stop band attenuation in dB 20 Filter order: 26 Filter coefficients are: b 0.0301 -0.0393 -0.0055 0.0172 -0.0285 -0.0117 0.0390 -0.0291 -0.0330 0.0759 -0.0308 -0.1191 0.2891 0.6352 0.2891 -0.1191 -0.0308 0.0759 -0.0330 -0.0291 0.0390 -0.0117 -0.0285 0.0172 -0.0055 -0.0393 0.0301


Page 115

08.508


2. FIR high pass filter using optimum equiripple FIR filter design using optimum equiripple algorithm Enter filter type (LPF, HPF, BPF or BSF): 'HPF' Enter pass band frequency fp in Hz 5000 Enter stop band frequency fs in Hz 4500 Enter sampling frequency Fs in samples/sec 15000 Enter pass band ripple in dB 1 Enter stop band attenuation in dB 20 Filter order: 26 Filter coefficients are: b 0.0050 0.0563 -0.0003 -0.0242 0.0253 0.0063 -0.0404 0.0305 0.0288 -0.0771 0.0341 0.1160 -0.2911 0.3687 -0.2911 0.1160 0.0341 -0.0771 0.0288 0.0305 -0.0404 0.0063 0.0253 -0.0242 -0.0003 0.0563 0.0050

3. FIR band pass filter using optimum equiripple FIR filter design using optimum equiripple algorithm Enter filter type (LPF, HPF, BPF or BSF): 'BPF' Enter stop band frequencies fs1 and fs2 (fs2 >fs1) in Hz [4000 7000] Enter pass band frequencies fp1 and fs2 (fp1 >fs1 and fp2

Page 116

08.508


4. FIR band stop filter using optimum equiripple FIR filter design using optimum equiripple algorithm Enter filter type (LPF, HPF, BPF or BSF): 'BSF' Enter stop band frequencies fs1 and fs2 (fs2 >fs1) in Hz [3000 5000] Enter pass band frequencies fp1 and fs2 (fp1 < fs1 and fp2 > fs2) in Hz [2500 5500] Enter sampling frequency Fs in samples/sec 15000 Enter pass band ripple in dB 1 Enter stop band attenuation in dB 30 Filter order: 36 Filter coefficients are: b 0.0011 -0.0079 -0.0003 -0.0461 0.0027 0.0206 -0.0010 0.0189 -0.0221 -0.0494 0.0423 0.0302 -0.0033 0.0286 -0.1204 -0.0645 0.2699 0.0336 0.6622 0.0336 0.2699 -0.0645 -0.1204 0.0286 -0.0033 0.0302 0.0423 -0.0494 -0.0221 0.0189 -0.0010 0.0206 0.0027 -0.0461 -0.0003 -0.0079 0.0011


Page 117

08.508 Expt. No. 9

Digital Signal Processing Lab STUDY OF SAMPLING RATE CONVERSION BY A RATIONAL FORM USING MATLAB

Objectives: To study sampling rate conversion by a rational form using MATLAB Equipments required: PC installed with MATLAB 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 in serting 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 of L • 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:

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. Department of ECE, VKCET

Page 118

08.508


Interpolator : To increase the sampling rate by an integer factor L. If x a(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 x i(n), except those that are at integer multiples of 2π. This is accomplished by filtering x i(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.

SRC by rational factor: SRC by L/M requires performing an interpolation to a sampling rate which is divisible by both L and M. The final output is then achieved by decimating by a factor of M. The need for a non-integer sampling rate conversion appears when the two systems operating at different sampling rates have to be connected, or when there is a need to convert the sampling rate of the recorded data into another sampling rate for further processing or reproduction. Such applications are very common in telecommunications, digital audio, multimedia and others. An example is transferring data from compact disc (CD) system at a rate of 44.1 kHz to a digital audio tape at 48 kHz. This can be achieved by increasing the data rate of the CD by a factor of 48/44.1, a non-integer.


Page 119

08.508


Illustration for sampling rate converter is:

If M>L, the resulting operation is a decimation process by a non-integer, and when M

Page 120

08.508


Programs: %Sampling rate conversion by rational factor clc; clear all; close all; disp('Sampling rate conversion by rational factor '); f=input('Enter frequency of the signal '); Fs=input('Enter sampling rate in samples/sec '); N=input('Enter number of samples N = '); L=input('Enter up-sampling factor L = '); M=input('Enter down-sampling factor M = '); %Time index Ts=1/Fs; n=0:Ts:(N-1)*Ts; %Amplitude A=2; %Input sequence x=A*sin(2*pi*f*n); %Up-sampling xi=interp(x,L); %Down-sampling xd=decimate(x,M); %Re-sampling y=resample(x,L,M); %Plotting sequences %Input sequence n1=0:N-1; subplot(4,1,1); stem(n1,x); xlabel('Time index n'); ylabel('x(n)'); title(['Input signal with sampling rate Fs = ',num2str(Fs), 'Hz']); %Up-sampled sequence ni=0:length(xi)-1; subplot(4,1,2); stem(ni,xi); xlabel('Time index n'); ylabel('xi(n)'); Fsu=Fs*L; title(['Up sampled signal with sampling rate Fs = ',num2str(Fsu),' Hz']); %Down-sampled sequence nd=0:length(xd)-1; subplot(4,1,3); stem(nd,xd); xlabel('Time index n'); ylabel('xd(n)'); Department of ECE, VKCET

Page 121

08.508


Fsd=Fs/M; title(['Down sampled signal with sampling rate Fs = ',num2str(Fsd),' Hz']); %Re-sampled sequence ny=0:length(y)-1; subplot(4,1,4); stem(ny,y); xlabel('Time index n'); ylabel('y(n)'); Fssc=Fs*L/M; title(['Sampling rate converted signal with sampling rate L/M = ',num2str(Fssc),' Hz']); Results: 1. Sampling rate conversion by rational factor Enter frequency of the signal 200 Enter sampling rate in samples/sec 4000 Enter number of samples N = 100 Enter up-sampling factor L = 1 Enter down-sampling factor M = 3

2. Sampling rate conversion by rational factor Enter frequency of the signal 200 Enter sampling rate in samples/sec 4000 Enter number of samples N = 100 Enter up-sampling factor L = 3 Enter down-sampling factor M = 1


Page 122

08.508





Page 123

08.508

Digital Signal Processing Lab BIBLIOGRAPHY

Text Books: 1. Vinay K. Ingle & John Proakis - Digital Signal Processing A MATLAB based Approach 2. Sanjit K. Mitra - Digital Signal Processing – Computer Based Approach 3. Emmanuel C. Ifeachor & Barrie W, Jervis – Digital Signal Processing A practical approach

References:

1. John G. Proakis & Dimitris K Manolakis - Digital Signal Processing 2. Alan V. Oppenheim Ronald W. Schafer - Discrete-Time Signal Processing 3. Steven W. Smith - The Scientist & Engineer's Guide to Digital Signal Processing 4. Edmund Lai – Practical Digital Signal Processing 5. E.S. Gopi - Algorithm Collections for Digital Signal Processing Applications using Matlab 6. Catherine Forbes, Merran Evans, Nicholas Hastings, and Brian Peacock Statistical Distributions


Page 124

08.508 DSP Lab Manual Part-B

Recommend Documents