Dsp First Lab Manual 1464v

Lab Manual

CSE‐205 Introduction to DSP       

         

 

Prepared by Fariha Atta Department of Computer Systems Engineering N-W.F.P University of Engineering and Technology Peshawar [email protected]

September 2008

         

 

Week-Wise Course Distribution

1st Week 

2nd Week 

3rd Week 

4th Week 

5th Week 

6th Week 

Introduction  Signals.           Different Types of Signals  Systems           Different Types of Systems  Mathematical representation of signals  Mathematical representation of systems  Examples of signals and systems  Continuous time and discrete time signals  Continuous time and discrete time systems  Sinusoids  Basic trigonometry  Sine and Cosine functions  Sinusoidal signals  Relation of frequency to period  Relation of phase shift to time shift  Sampling and plotting sinusoids  Complex exponentials and phasors    Review of complex numbers    Complex exponential signals    The rotating phasor interpretation    Inverse Euler formulas  Phasor addition    Addition of complex numbers    Phasor addition rule  Spectrum Representation  The spectrum of a sum of sinusoids    Graphical plot of spectrum  Beat notes    Multiplication of sinusoids    Beat note waveform    Amplitude modulation  Periodic waveforms  Fourier analysis  The square wave & The triangular wave  Non periodic signals  Time frequency spectrum  Frequency modulation: Chirp signals  Stepped Frequency 

7th Week 

8th Week 

  9th Week   

10th   Week 

11th Week 

12th   Week 

Sampling and Aliasing  Sampling  Sampling sinusoidal signals  The sampling theorem  Aliasing              Folding    Spectrum view of sampling  Over‐sampling  Aliasing due to under sampling  Folding due to under sampling  Discrete‐to‐continuous conversion    Course Review    Mid Term Exam  FIR Filters  Discrete time systems  The running average filter  The general FIR filter    An illustration of FIR filter  The unit impulse response  Unit impulse sequence  Unit impulse response sequence  The unit delay system  Convolution and FIR filters    Computing the output of a convolution  Implementation of FIR Filters              Block Diagrams  Linear Time‐Invariant systems  Time invariance  Linearity  The FIR case  Convolution and LTI systems  Some properties of LTI systems    Convolution as an operator    Commutative property of Convolution    Associative property of Convolution  Cascaded LTI systems   

13th Week 

14th Week 

15th   Week 

16th   Week 

Frequency response of FIR filters  Sinusoidal response of FIR filters  Superposition and Frequency Response  Properties of the Frequency Response  Graphical representation of the frequency response              Delay System              First Difference System              A simple low  filter  Cascaded LTI Systems  Running average filtering              Plotting the Frequency Response    Z‐Transforms  Definition of the Z‐Transform  The Z‐Transform and linear systems  Properties of the Z‐Transform  A general Z‐Transform formula  The Z‐Transform as an operator  Unit delay operator  Operator notation  Operator notation in block diagram  Convolution and the Z‐Transform            Cascading Systems            Factorizing Z‐Polynomials            Deconvolution  Relationship between the Z‐Domain and w‐Domain            The Z‐Plane and the Unit Circle            The Zeros and Poles of H(z)             Significance of the Zeros of H(z)  Inverse Z‐Transform  IIR Filters  The general IIR difference equation  Time domain response  Linear and time invariance of IIR filters  Impulse response of a first order IIR system  System Function of an IIR Filter           The  General First order Case  Poles and Zeros  Poles or zeros at the origin or infinity  Pole locations and stability  Frequency response using Matlab 

17th   Week 

Spectrum analysis  Review of the frequency spectrum  Spectrum analysis of Periodic signals    Periodic signals 

Spectrum of a periodic signal Course Review (1 Week)   18th   Week 

Final Exam  

  Objectives of Lab Course of DSP 1st  In  the  lab,  the  students  will  acquire  hands‐on  experience  with  programming  in  MATLAB.  MATLAB  will  enable  them  to  study  and  understand  the  basics  of  Digital  Signal  Processing  as  well  as  validate  the  theory with real‐world examples. The labs will cover complex exponential signal, synthesize complicated  sinusoidal  waveforms,  response  of  FIR  filters,  sampling,  along  with  several  interesting  digital  signal‐ processing  (DSP)  applications.  For  the  Lab  part  grading  will  largely  depend  on  hands‐on  proficiency  of  the students in DSP related portion of the MATLAB 7.0. 

 

Lab # 1

OBJECTIVES OF THE LAB

---------------------------------------------------------------------Matlab will be used extensively in all the succeeding labs. The goal of this first lab is to gain familiarity with Matlab and build some basic skills in the Matlab language. Some specific topics covered in this lab are: • Introduction to Matlab • Matlab Environment • Matlab Help • Variable arithmetic • Built in Mathematical Functions • Input and display • Timing functions • Introduction to M-files

----------------------------------------------------------------------

1.1 WHAT IS MATLAB? MATLAB is a commercial "MATrix LABoratory" package, by MathWorks, which operates as an interactive programming environment with graphical output. The MATLAB programming language is exceptionally straightforward since almost every data object is assumed to be an array. Hence, for some areas of engineering MATLAB is displacing popular programming languages, due to its interactive interface, reliable algorithmic foundation, fully extensible environment, and computational speed.

1.2 ENTERING AND RUNNING MATLAB Double click on the MATLAB icon to launch and a command window will appear with the prompt: >> You are now in MATLAB. From this point on, individual MATLAB commands may be given at the program prompt. They will be processed when you hit the <ENTER> key. The following figure shows the screenshot of matlab.

1.3 LEAVING MATLAB A MATLAB session may be terminated by simply typing >> quit or by typing >> exit at the MATLAB prompt.

1.4 MATLAB HELP Online help is available from the MATLAB prompt, both generally (listing all available commands). >> help [a long list of help topics follows] and for specific commands: >> help command_name

If you want to search for all the commands related to some particular functionality, use the keyword lookfor followed by a keyword that explains the functionality. >>lookfor convolution will return a number of commands that perform convolution related tasks.

1.5 VARIABLES MATLAB has built-in variables like pi, eps, and ans. You can learn their values from the MATLAB interpreter. >> eps eps = 2.2204e-16 >> pi ans = 3.1416

1.5.1

Variable Assignment

The equality sign is used to assign values to variables: >> x = 3 x= 3 >> y = x^2 y= 9 Variables in MATLAB are case sensitive. Hence, the variables "x" and "y" are distinct from "X" and "Y" (at this point, the latter are in fact, undefined). Output can be suppressed by appending a semicolon to the command lines. >> x = 3; >> y = x^2; >> y

y= 9

1.5.2

Active Variables

At any time you want to know the active variables you can use who: >> who Your variables are: ans x y

1.5.3

Removing a Variable

To remove a variable, try this: >> clear x To remove all the variables from workspace, use clear >> clear

1.5.4

Saving and Restoring Variables

To save the value of the variable "x" to a plain text file named "x.value" use >> save x.value x -ascii To save all variables in a file named mysession.mat, in reloadable format, use >> save mysession To restore the session, use >> load mysession

1.6 VARIABLE ARITHMETIC MATLAB uses some fairly standard notation. More than one command may be entered on a single line, if they are separated by commas. >> 2+3; >> 3*4, 4^2;

Powers are performed before division and multiplication, which are done before subtraction and addition. For example >> 2+3*4^2; generates ans = 50. That is: 2+3*4^2 ==> 2 + 3*4^2 <== exponent has the highest precedence ==> 2 + 3*16 <== then multiplication operator ==> 2 + 48 <== then addition operator ==> 50

1.6.1

Double Precision Arithmetic

All arithmetic is done to double precision, which for 32-bit machines means to about 16 decimal digits of accuracy. Normally the results will be displayed in a shorter form. >> a = sqrt(2) a= 1.4142 >> format long, b=sqrt(2) b= 1.41421356237310 >> format short

1.6.2

Command-Line Editing

The arrow keys allow "command-line editing," which cuts down on the amount of typing required, and allows easy error correction. Press the "up" arrow, and add "/2." What will this produce? >> 2+3*4^2/2 Parentheses may be used to group , or to make them more readable. For example: >> (2 + 3*4^2)/2

generates ans = 25.

1.6.3

Built-In Mathematical Functions

MATLAB has a platter of built-in functions for mathematical and scientific computations. Here is a summary of relevant functions.

Function Meaning Example ====================================================== sin sine sin(pi) = 0.0 cos cosine cos(pi) = 1.0 tan tangent tan(pi/4) = 1.0 asin arcsine asin(pi/2)= 1.0 acos arccosine acos(pi/2)= 0.0 atan arctangent atan(pi/4)= 1.0 exp exponential exp(1.0) = 2.7183 log natural logarithm log(2.7183) = 1.0 log10 logarithm base 10 log10(100.0) = 2.0 ======================================================

The arguments to trigonometric functions are given in radians.

Example: Let's that sin(x)^2 + cos(x)^2 = 1.0 for arbitrary x. The MATLAB code is: >> x = pi/3; >> sin(x)^2 + cos(x)^2 - 1.0 ans = 0

1.7 TIMING COMMANDS Timing functions may be required to determine the time taken by a command to execute or an operation to complete. Several commands are available to accomplish it:

1.7.1

Clock

CLOCK returns Current date and time as date vector.

CLOCK returns a six

element date vector vector containing the current time and date in decimal form: CLOCK = [year month day hour minute seconds] The first five elements are integers. The second’s element is accurate to several digits beyond the decimal point. FIX(CLOCK) rounds to integer display format.

1.7.2

Etime

ETIME Elapsed time. ETIME(T1,T0) returns the time in seconds that has elapsed between vectors T1 and T0. The two vectors must be six elements long, in the format returned by CLOCK: T = [Year Month Day Hour Minute Second] Time differences over many orders of magnitude are computed accurately. The result can be thousands of seconds if T1 and T0 differ in their first five components or small fractions of seconds if the first five components are equal. t0 = clock; operation etime(clock,t0)

1.7.3

Tic Toc

TIC Start a stopwatch timer. The sequence of commands TIC, operation, TOC Prints the number of seconds required for the operation.

1.8 INPUT & DISPLAY 1.8.1

INPUT

INPUT prompts for input. R = INPUT('How many apples')

gives the the prompt in the text string and then waits for input from the keyboard. The input can be any MATLAB expression, which is evaluated, using the variables in the current workspace, and the result returned in R. If the presses the return key without

entering anything, INPUT returns an empty matrix.

Example: Entering a single variable >> x=input('Enter a variable: ') Enter a variable: 4 x= 4 >> x=input('Enter a vector: ')

Example: Entering a vector A vector is entered by specifying [] and elements are inserted inside these brackets, separated by space. Enter a vector: [3 4 1] x= 3

4

1

Example: A \n entered after the string results in starting a new line. >> x=input('Enter a value\n') Enter a value 5 x= 5

1.8.2

DISP

DISP Display array.

DISP(X) displays the array, without printing the array name. In all other ways it's the same as leaving the semicolon off an expression except that empty arrays don't display. DISP(‘string’) is another variation of the same function that is used to display a string on the command prompt. Example: >> disp('I am using MATLAB 7.0') I am using MATLAB 7.0

1.9 M-Files Typing errors are time-consuming to fix if you are working in the command window because you need to retype all or part of the program. Even if you do not make any mistakes, all of your work may be lost if you inadvertently quit MATLAB. To preserve large sets of commands, you can store them in a special type of file called an M-file. MATLAB s two types of M-files: script and function M-files. To hold a large collection of commands, we use a script M-file. The function M-file is discussed in coming lab. The script file has a '.m' extension and is referred to as an M-file (for example, myfile.m myfuncion.m, etc.). The commands in the script file can then be executed by typing the file name without its extension in the command window. Commands in a script utilize and modify the contents of the current workspace. It is possible to embed comments in a script file. To make a script M-file, you need to open a file using the built-in MATLAB editor. There are two ways to accomplish it: 1. From file menu, click NEW 2. Type edit on command line A new window appears like one shown in the figure below.

When you are finished with typing in this new window, click File->Save to save this file. The extension of this file be .m. In order to execute this program, 1. Write the name of file on command window (excluding the .m) or 2. Click Debug->Run

---------------------------TASK 1---------------------------Create an m-file and write a program for calculating area of a circle. Try out several other programs of similar computation.

---------------------------TASK 2----------------------------

Create an m-file to get 10 numbers from and generate the square of those numbers.

Lab # 2

OBJECTIVES OF THE LAB

---------------------------------------------------------------------In this lab, we will cover the following topics: • • • • •

Built in Matrix Functions Indexing Matrices Sub Matrices Matrix element level operations Round Floating Point numbers to Integers

----------------------------------------------------------------------

1.1 MATRICES MATLAB works with essentially only one kind of object, a rectangular numerical matrix possibly, with complex entries. Every MATLAB variable refers to a matrix [a number is a 1 by 1 matrix]. In some situations, 1-by-1 matrices are interpreted as scalars, and matrices with only one row or one column are interpreted as vectors. A matrix is a rectangular array of numbers. For example:

⎡3 ⎢1 ⎢ ⎢2 ⎢ ⎣1

6 4 8 4

9 8 7 2

2⎤ 5⎥⎥ 5⎥ ⎥ 3⎦

defines a matrix with 3 rows, 4 columns, and 12 elements.

Example: consider the following three equations: 3 * x1 - 1 * x2 + 0 * x3 = 1 -1 * x1 + 4 * x2 - 2 * x3 = 5 0 * x1 - 2 * x2 + 10 * x3 = 26 This family of equations can be written in the form A.X = B, where [ 3 -1 0 ] A=

[-1 4 -2 ], [0 -2 10 ]

[ x1 ] X=

[ x2 ], and [ x3 ]

[1] B=

[5] [ 26 ]

Depending on the specific values of coefficients in matrices A and B, there may be: (a) no solutions to A.X = B, (b) a unique solution to A.X = B, or (c) an infinite number of solutions to A.X = B. In this particular case, however, the solution matrix [1] X=

[2] [3]

makes the right-hand side of the matrix equations (i.e., A.X) equal the left-hand side of the matrix equations (i.e., matrix B).

1.1.1

Defining Matrices In Matlab

MATLAB is designed to make definition of matrices and matrix manipulation as simple as possible. Matrices can be introduced into MATLAB in several different ways: For example, either of the statements >> A = [1 2 3; 4 5 6; 7 8 9]; and >> A = [ 1 2 3 456 789] creates the obvious 3-by-3 matrix and assigns it to a variable A. Note that: •

The elements within a row of a matrix may be separated by commas as well as a blank.

•

The elements of a matrix being entered are enclosed by brackets;

•

A matrix is entered in "row-major order" [i.e. all of the first row, then all of the second row, etc];

•

Rows are separated by a semicolon [or a newline], and the elements of the row may be separated by either a comma or a space. [Caution: Watch out for extra spaces!]

The matrix element located in the i-th row and j-th column of a is referred to in the usual way: >> A(1,2), A(2,3) ans = 2 ans = 6 It's very easy to modify matrices: >> A(2,3) = 10;

1.1.2

Building Matrices from a Block

Large matrices can be assembled from smaller matrix blocks. For example, with matrix A in hand, we can enter the following commands: >> C = [A; 10 11 12]; <== generates a (4x3) matrix >> [A; A; A]; <== generates a (9x3) matrix >> [A, A, A]; <== generates a (3x9) matrix As with variables, use of a semicolon with matrices suppresses output. This feature can be especially useful when large matrices are being generated.

1.1.3

Built-in matrix functions

MATLAB has many types of matrices which are built into the system e.g. Function Description =============================================== diag returns diagonal M.E. as vector eye identity matrix hilb Hilbert matrix magic magic square ones matrix of ones rand randomly generated matrix triu upper triangular part of a matrix tril lower triangular part of a matrix zeros matrix of zeros ===============================================

Here are some examples: i. Matrices of Random Entries: A 3 by 3 matrix with random entries is produced by typing >> rand(3) ans = 0.0470

0.9347

0.8310

0.6789

0.3835

0.0346

0.6793

0.5194

0.0535

General m-by-n matrices of random entries are generated with >> rand(m,n);

ii.

Magic Squares: A magic square is a square matrix which has equal sums along all its rows and columns. For example: >> magic(4) ans = 16

2

3

13

5

11

10

8

9

7

6

12

4

14

15

1

The elements of each row and column sum to 34. iii.

Matrices of Ones: The functions eye (m,n) produces an m-by-n matrix of ones. eye (n) produces an n-by-n matrix of ones.

iv.

Matrices of Zeros: The commands zeros (m,n) produces an m-by-n matrix of zeros. zeros (n) produces an n-by-n one; If A is a matrix, then zeros (A) produces a matrix of zeros of the same size as A.

v.

Diagonal Matrices: If x is a vector, diag(x) is the diagonal matrix with x down the diagonal. If A is a square matrix, then diag(A) is a vector consisting of the diagonal of A. What is diag(diag(A))? Try it.

1.2 MATRIX OPERATIONS The following matrix operations are available in MATLAB: Operator Description Operator Description ============================================================ + addition ' transpose subtraction \ left division * multiplication / right division ^ power ============================================================

These matrix operations apply, of course, to scalars (1-by-1 matrices) as well. If the sizes of the matrices are incompatible for the matrix operation, an error message will result, except in the case of scalar-matrix operations (for addition, subtraction, and division as well as for multiplication) in which case each entry of the matrix is operated on by the scalar.

1.2.1

Matrix Transpose

The transpose of a matrix is the result of interchanging rows and columns. MATLAB denotes the [conjugate] transpose by following the matrix with the singlequote [apostrophe]. For example: >> A' ans = 1

4

7

2

5

8

3

6

9

>> B = [1+i

2 + 2*i

3 - 3*i];

>> B = B' B=

1.2.2

1.0000

- 1.0000i

2.0000

- 2.0000i

3.0000

+ 3.0000i

Matrix Addition/Subtraction

Let matrix "A" have m rows and n columns, and matrix "B" have p rows and q columns. The matrix sum "A + B" is defined only when m equals p and n equals q, the result is a n-by-m matrix having the element-by-element sum of components in A and B. For example: >> A = [ 2 3; 4 5.0; 6 7]; >> B = [ 1 -2; 3 6.5; 10 -45]; >> A+B

ans =

1.2.3

3.0000

1.0000

7.0000

11.5000

16.0000

-38.0000

Matrix Multiplication

Matrix multiplication requires that the sizes match. If they don't, an error message is generated. >> A*B, B*A; <== results into error as inner matrix dimensions doesn’t agrees >> B'*A; >> A*A', A'*A; >> B'*B, B*B'; Scalars multiply matrices as expected, and matrices may be added in the usual way (both are done "element by element”): >> 2*A, A/4; >> A + [b,b,b]; <== results into error as inner matrix dimensions doesn’t agrees

Example: We can use matrix multiplication to check the "magic" property of magic squares. >> A = magic(5); >> b = ones(5,1); >> A*b;

<== (5x1) matrix containing row sums.

>> v = ones(1,5); >> v*A;

1.2.4

<== (1x5) matrix containing column sums.

Matrix Functions "any" and "all”

There is a function to determine if a matrix has at least one nonzero entry, any, as well as a function to determine if all the entries are nonzero, all. >> A = zeros(1,4) >> any(A) >> D = ones(1,4) >> any(D) >> all(A)

1.2.5

Returning more than One Value

Some MATLAB functions can return more than one value. In the case of max the interpreter returns the maximum value and also the column index where the maximum value occurs. Similarly, min function returns the minimum value along with the column index where the minimum value occurs. >> B = magic(4); >> [m, i] = max(B) >> [m, i] = min(B)

1.2.6

Size of Matrix

Size of a matrix can be calculate by using function ‘size ‘. >> x = [1 2 3 ;1 2 3]; >> s = size(x) s= 2

1.2.7

3

Length of Array

Length of an array can be found using function ‘length’. >> n = [-3:1:3]; >> l = length(n) l= 7

1.2.8

Finding an element in a matrix

This function can be used to find index of any particular value. Say given array is >> x= [0 2 4 6 8]; To find the indices of all values that are greater than 4, following is used >> y = find(x>4) y= 4

5

---------------------------TASK 1-----------------------------

Write a program to generate a new matrix from the matrix given below such that each column in the new matrix except the first one is the result of subtraction of that column from the previous one i.e. 2nd new column is the result of subtraction of 2nd column and 1st column and so on. Copy the first column as it is in the new matrix. ⎡3 ⎢1 ⎢ ⎢2 ⎢ ⎣1

6 4 8 4

9 8 7 2

2⎤ 5⎥⎥ 5⎥ ⎥ 3⎦

---------------------------TASK 2-----------------------------

Generate two 10000 sampled random discrete time signals (1 dimensional) using rand() function i.e. rand(1,10000). Write a program to add the two signals together using simple vector addition. Determine the time required for addition using tic, toc pair or etime function.

1.3 SUB-MATRICES A note about Colon Notation A central part of the MATLAB language syntax is the "colon operator," which produces a list. For example: >> -3:3 ans = -3 -2 -1 0 1 2 3 The default increment is by 1, but that can be changed. For example: >> x = -3 : .3 : 3 x= Columns 1 through 7 -3.0000 -2.7000 -2.4000 -2.1000 -1.8000 -1.5000 -1.2000 Columns 8 through 14 -0.9000 -0.6000 -0.3000 0 0.3000 0.6000 0.9000 Columns 15 through 21 1.2000 1.5000 1.8000 2.1000 2.4000 2.7000 3.0000

This can be read: "x is the name of the list, which begins at -3, and whose entries increase by .3, until 3 is sured." You may think of x as a list, a vector, or a matrix, whichever you like. In our third example, the following statements generate a table of sines. >> x = [0.0:0.1:2.0]'; >> y = sin(x); >> [x y] Try it. Note that since sin operates entry-wise, it produces a vector y from the vector x. The colon notation can also be combined with the earlier method of constructing matrices. >> a = [1:6 ; 2:7 ; 4:9]

---------------------------TASK 3-----------------------Generate the sequence -99, -96, -93, . . ., -3, 0, 3, 6 , . . . .,93, 96, 99.

Colon notation can be used to generate vectors. A very common use of the colon notation is to extract rows, or columns, as a sort of "wild-card" operator which produces a default list. For example, A(1:4,3)

is the column vector consisting of the first four entries of the third column of A .

A(:,3)

is the third column of A. A colon by itself denotes an entire row or column.

A(1:4,:)

is the first four rows of A.

Arbitrary integral vectors can be used as subscripts. The statement A(:,[2 4])

contains as columns, columns 2 and 4 of matrix A.

This subscripting scheme can be used on both sides of an assignment statement: A(:,[2 4 5]) = B(:,1:3) replaces columns 2,4,5 of matrix A with the first three columns of matrix B. Note that the "entire" altered matrix A is printed and assigned. Try it.

---------------------------TASK 4------------------------Create two matrices i.e. A consisting of 1 through 6 & 12 through 7, while B consisting of 6 through 1 & 7 through 12. Perform the following operations: A+B,

A-B, A.*B, A./B, A.^2, 1./A, A/2, A+1. Take matrices of your choice and perform the above mentioned operations on them.

---------------------------TASK 5-------------------------

MATLAB has functions to round floating point numbers to integers. These are round, fix, ceil, and floor. Test how these functions work. Determine the output of the following: >> f = [-.5 .1 .5] >> round(f) >> fix(f) >> ceil(f) >> floor(f) >> sum(f) >> prod(f)

Lab # 3

OBJECTIVES OF THE LAB

---------------------------------------------------------------------In this lab, we will get an understanding of the following topics: • • • • • •

Making Functions Control Structures Relational Constructs Logical Constructs Branching Constructs Looping constructs

----------------------------------------------------------------------

1.1 MAKING FUNCTIONS A function can be created by the following syntax: function [output1, output2, ...] = cmd_name(input1,input2,...) A function is a reusable portion of code that can be called from program to accomplish some specified functionality. A function takes some input arguments and returns some output. To create a function that adds two numbers and stores the result in a third variable, type the following code in an m-file: function add x=3; y=5; z=x+y Save the file by the name of add (in work folder, which is chosen by default), go back to the command window and write >> add z= 8 You see that the sum z is displayed in the command window. Now go back to the editor/debugger and modify the program as follows function addv(x,y) z=x+y Save the above program with a new name addv, go back to the command window and type the following >> addv(3,5) z= 8 >> addv(5,5) z= 10 We have actually created a function of our own and called it in the main program and gave values to the variables (x,y).

Once again go back to the editor/debugger and modify the program as follows function adv(x,y) %------------------------------------------------% This function takes two values as input, % finds its sum, & displays the result. % inputs: x & y % output: z % Example: addv(3,6) % Result: z=9 %-------------------------------------------------z=x+y Save the program with the same name adv, go back to command window, type the following >> help addv ------------------------------------------------This function takes two values as input, finds its sum, & displays the result. inputs: x & y output: z Example: addv(3,6) Result: z=9 -------------------------------------------------SCRIPT VS FUNCTION •

A script is simply a collection of Matlab commands in an m-file. Upon typing the name of the file (without the extension), those commands are executed as if they had been entered at the keyboard. Functions are used to create -defined matlab commands.

•

A script can have any name. A function file is stored with the name specified after keyword function.

•

The commands in the script can refer to the variables already defined in Matlab, which are said to be in the global workspace. When a function is invoked, Matlab creates a local workspace. The commands in the function cannot refer to variables from the global (interactive) workspace

unless they are ed as inputs. By the same token, variables created as the function executes are erased when the execution of the function ends, unless they are ed back as outputs.

---------------------------TASK 1---------------------------Construct a function in M-file by the name of greater(x,y), which will take two inputs from the , finds the value that is greater among the two and then displays it.

1.2 CONTROL STRUCTURES Control-of-flow in MATLAB programs is achieved with logical/relational constructs, branching constructs, and a variety of looping constructs.

1.2.1

Relational and logical constructs

The relational operators in MATLAB are Operator Description =================================== < less than > greater than <= less than or equal >= greater than or equal == equal ~= not equal ===================================

Note that ``='' is used in an assignment statement while ``=='' is used in a relation. Relations may be connected or quantified by the logical operators Operator Description =================================== & and | or ~ not ===================================

When applied to scalars, a relation is actually the scalar 1 or 0 depending on whether the relation is true or false (indeed, throughout this section you should think of 1 as true and 0 as false). For example >> 3 < 5

ans = 1 >> a = 3 == 5 a= 0 When logical operands are applied to matrices of the same size, a relation is a matrix of 0's and 1's giving the value of the relation between corresponding entries. For example: >> A = [ 1 2; 3 4 ]; >> B = [ 6 7; 8 9 ]; >> A == B ans = 0

0

0

0

>> A < B ans = 1

1

1

1

To see how the other logical operators work, you should also try >> ~A >> A&B >> A & ~B >> A | B >> A | ~A

1.2.2

Branching constructs

MATLAB provides a number of language constructs for branching a program's control of flow.

i.

if-end Construct : The most basic construct is if <program> end Here the condition is a logical expression that will evaluate to either true or false (i.e., with values 1 or 0). When the logical expression evaluates to 0, the program control moves on to the next program construction. You should keep in mind that MATLAB regards A==B and A<=B as functions with values 0 or 1.

Example: >> a = 1; >> b = 2; >> if a < b c = 3; end; >> c c= 3 ii.

If-else-end Construct: Frequently, this construction is elaborated with if <program1> else <program2> end In this case if condition is 0, then program2 is executed.

iii.

If-elseif-end Construct: Another variation is if <program1> elseif <program2> end

Now if condition1 is not 0, then program1 is executed, if condition1 is 0 and if condition2 is not 0, then program2 is executed, and otherwise control is ed on to the next construction.

----------------------------TASK 2--------------------------Find for integer 0 < a ≤ 10, the values of C, defined as follows: C= 5ab, C= ab,

0 <=5 5 <=10

where b = 15.

----------------------------TASK 3--------------------------For the values of integer a going from 1 to 10, using separately the methods of if syntax and the Boolean alternative expressions, find the values of C if: C= 2a, C= a + 5, C= a, 1.2.3 i.

a<3 3<=a<7 a rel="nofollow">7

Looping constructs For Loops : A for loop is a construction of the form for i= 1 : n <program> end

The program will repeat <program> once for each index value i = 1, 2, .... n. Here are some examples of MATLAB's for loop capabilities:

Example: The basic for loop >> for i = 1 : 5, c = 2*i end c= 2

..... lines of output removed ... c= 10 computes and prints "c = 2*i" for i = 1, 2, ... 5.

Example: For looping constructs may be nested. Here is an example of creating matrices contents inside a nested for loop: >> for i=1:10 for j=1:10 A(i,j) = i/j; end end There are actually two loops here, with one nested inside the other; they define A(1,1), A(1,2), A(1,3) ... A(1,10), A(2,1), ... A(10,10) in that order.

Example: MATLAB will allow you to put any vector in place of the vector 1:n in this construction. Thus the construction >> for i = [2,4,5,6,10] <program> end is perfectly legitimate. In this case program will execute 5 times and the values for the variable i during execution are successively, 2,4,5,6,10.

----------------------------TASK 4--------------------------Generate the square of the first ten integers.

----------------------------TASK 5---------------------------Add the following two matrices using for loop.

⎡5 12 3⎤ ⎢ 9 6 5⎥ ⎥ ⎢ ⎢⎣2 2 1⎥⎦

ii.

⎡ 2 1 9⎤ ⎢10 5 6⎥ ⎥ ⎢ ⎢⎣ 3 4 2⎥⎦

While Loops

A while loop is a construction of the form while <program> end where condition is a MATLAB function, as with the branching construction. The program will execute successively as long as the value of condition is not 0. While loops carry an implicit danger in that there is no guarantee in general that you will exit a while loop. Here is a sample program using a while loop. function l=twolog(n) % l=twolog(n). l is the floor of the base 2 % logarithm of n. l=0; m=2; while m<=n l=l+1; m=2*m; end

----------------------------TASK 6---------------------------Create an m-file that inputs a number from and then finds out the factorial of that number.

----------------------------TASK 7----------------------------

Create an m-file that takes two vectors from . Make sure that the second vector taken is of the same size as the first vector (Hint: use while loop). In a while loop, generate a third vector that contains the sum of the squares of corresponding entries of both the vectors.

Lab # 4

OBJECTIVES OF THE LAB

---------------------------------------------------------------------This lab will help you grasp the following concepts: • • • • •

Discrete Signal representation in Matlab Matlab Graphics Two Dimensional Plots Plot and subplot Different Plotting Functions Used in Matlab

----------------------------------------------------------------------

4.1 DISCRETE-TIME SIGNAL REPRESENTATION IN MATLAB In MATLAB, finite-duration sequence (or discrete time signal) is represented by row vector of appropriate values. Such representation does not have any information about sample position n. Therefore, for correct representation, two vectors are required, one for x and other for n.

Consider the following finite duration sequence & its

implementation: x(n) = { 1 -1 0 2 1 4 6 } ↑ >> n = [-3:1:3] n= -3

-2

-1

0

1

2

3

2

1

4

6

>> x = [1 -1 0 2 1 4 6] x= 1

-1

0

NOTE # 01: When the sequence begins at n=0, x-vector representation alone is enough. NOTE # 02: An arbitrary infinite-sequence can’t be represented in MATLAB due to limited memory.

-----------------------------TASK 1---------------------------Given the signals: X1[n] = [2 5 8 4 3] X2[n] = [4 3 2] a) Write a Matlab program that adds these two signals. Use vector addition and multiplication. Apply if-else construct, where condition in if-part checks the relative lengths of two vectors & performs the desired operations, otherwise in else-part it asks about two choices 1: exit from the program, 2: add redundant samples (equal to the difference of vector lengths) in the small vector, thereby creating new vector x2_mod. Use x2_mod to perform vector addition and multiplication. To implement this, use switch construct. b) Instead of using vector addition and multiplication, use for loop to add

and multiply the signals. Where for loop should run till the length of shortest sequence.

-----------------------------TASK 2---------------------------Amplitude scaling by a factor β causes each sample to get multiplied by β. Write a -defined function that has two input arguments: (i) a signal to be scaled and (ii) scaling factor β. The function should return the scaled output to the calling program. In the calling program, get the discrete time signal as well as the scaling factor from and then call the above-mentioned function.

-----------------------------TASK 3---------------------------X2[n] X1[n]

2

2

2

1

1

0

1

2

3

3

2 1

1

4

0 1

2

3

4

Write a Matlab program to compare the signals x1[n] and x2[n]. Determine the index where a sample of x1[n] has smaller amplitude as compared to the corresponding sample of x2[n]. Use for loop.

4.2 GRAPHICS Two- and three-dimensional MATLAB graphs can be given titles, have their axes labeled, and have text placed within the graph. The basic functions are: Function Description ============================================================================ plot(x,y) plots y vs x plot(x,y1,x,y2,x,y3) plots y1, y2 and y3 vs x on the same graph stem(x) plots x and draws a vertical line at each datapoint to the horizontal axis xlabel('x axis label') labels x axis ylabel('y axis label') labels y axis title ('title of plot') puts a title on the plot gtext('text') activates the use of the mouse to position a crosshair on the graph, at which point the 'text' will be placed when any key is pressed. zoom allows zoom IN/OUT using the mouse cursor grid draws a grid on the graph area print filename.ps saves the plot as a black and white postscript file Shg brings the current figure window forward. CLF clears current figure. ============================================================================

4.2.1

Two-dimensional plots

The plot command creates linear x-y plots; if x and y are vectors of the same length, the command plot(x,y) opens a graphics window and draws an x-y plot of the elements of x versus the elements of y. Example: Let's draw the graph of the sine function over the interval -4 to 4 with the following commands: >> x = -4:.01:4; y = sin(x); plot(x,y) >> grid; >> xlabel('x'); >> ylabel('sin(x)'); >> title('Graph of SINE function') The vector x is a partition of the domain with meshsize 0.01 while y is a vector giving the values of sine at the nodes of this partition (recall that sin operates entrywise). Following figure shows the result.

MULTIPLE PLOTS ON SAME FIGURE WINDOW Two ways to make multiple plots on a single graph are:

i.

Single plot command x = 0:.01:2*pi; y1=sin(x); y2=sin(2*x); y3=sin(4*x); plot(x,y1,x,y2,x,y3) xlabel('Time (sec)'); ylabel('Amplitude (A)');

ii.

Multiple plot commands

Another way is with hold. The command hold freezes the current graphics screen so that subsequent plots are superimposed on it. Entering hold again releases the ``hold.'' x = 0:.01:2*pi; y1=sin(x); y2=sin(2*x); y3=sin(4*x); plot(x,y1); hold on; plot(x,y2); plot(x,y3); xlabel('Time (sec)'); ylabel('Amplitude (A)');

OVERRIDING THE DEFAULT PLOT SETTINGS One can override the default linetypes and pointtypes. For example, the command sequence x = 0:.01:2*pi; y1=sin(x); y2=sin(2*x); y3=sin(4*x); plot(x,y1,'--',x,y2,':',x,y3,'+'); grid; title ('Dashed line and dotted line graph'); xlabel('Time (sec)'); ylabel('Amplitude (A)'); axis tight;

The line-type and mark-type are ============================================================= Linetypes : solid (-), dashed (--), dotted (:), dashdot (-.) Marktypes : point (.), plus (+), star (*), circle (o), x-mark (x) =============================================================

-------------------------TASK 4-------------------------Plot the two curves y1 = 2x + 3 and y2 = 4x + 3 on the same graph using different plot styles. AXES COMMANDS (MANUAL ZOOMING) MATLAB automatically adjusts the scale on a graph to accommodate the coordinates of the points being plotted. The axis scaling can be manually enforced by using the command axis([xmin xmax ymin ymax]). A signal can be zoomed out by specifying the axis coordinates by himself. Example: x = -5*pi:.01:5*pi; y1= sin(x); plot(x,y1,'r')

The plot is shown in the figure below.

In order to see only one cycle of this signal from 0 to 2π, the signal is zoomed using axis command. Here we have specified xmin and xmax as 0 and 2π respectively. x = -5*pi:0.01:5*pi; y1=sin(x); plot(x,y1, 'r') axis([0 2*pi -1 1]) The magnified plot is shown in the figure below.

Similarly the y-axis can be adjusted according to requirements. x = -5*pi:0.01:5*pi; y1=sin(x); plot(x,y1, 'r') axis([0 2*pi -2 2])

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('points in a plane')

SUBPLOT SUBPLOT Create axes in tiled positions. MATLAB graphics windows will contain one plot by default. The command subplot can be used to partition the screen so that up to four plots can be viewed simultaneously. A single figure can be divided into a number of plotting areas where different graphs can be plotted. This can be accomplished by using the command subplot(m, n, p) where m, n specifies the total number of rows

and columns respectively in the figure window and p specifies the specific cell to plot into. x = 0:1:10; y = x.^2; z = 10*x; Now type the following code figure subplot (1,2,1) plot(x,y) subplot (1,2,2) plot(x,z)

In the above case subplot(m,n,p) command was used, in our case subplot (1,2,1) and subplot (1,2,2). Here m=1 means that divide the figure into 1 row, n=2 means to divide the figure into 2 columns. This gives us a total of 2 subplots in one figure. Where p=1 means the window on the left (starting from row 1 and counting p=1 subplots to the right) and p=2 means the subplot on the right (starting from row 1 and counting p=2 subplots to the right).

Example: Performing operations on signals entered by clc clear all close all x = input('Enter the first discrete time signal\n'); len_x = length(x); y = input('Enter the second discrete time signal\n'); len_y = length(y); while(len_y~=len_x) disp('Error: Length of signals must match. Enter the 2nd signal again') y=input(''); len_y=length(y); end z = x+y; subplot(3,1,1); stem(x,'filled'); title('Signal 1'); xlabel('Sample number'); ylabel('Signal Amplitude'); subplot(3,1,2); stem(y,'filled'); title('Signal 2'); xlabel('Sample number'); ylabel('Signal Amplitude'); subplot(3,1,3); stem(z,'filled'); title('Resultant Signal'); xlabel('Sample number'); ylabel('Signal Amplitude');

output: Enter the first discrete time signal [3 5 1 0 2]

Enter the second discrete time signal [1 1 3 2 1]

------------------------------TASK 5---------------------------Make two separate functions for signal addition and multiplication. The functions should take the signals as input arguments and return the resultant signal. In the main program, get the signals from , call the functions for signal addition and multiplication, and plot the original signals as well as the resultant signals.

------------------------------TASK 6---------------------------Given the signals: X1[n] = 2δ[n] + 5δ[n-1] + 8δ[n-2] + 4δ[n-3] + 3δ[n-4] X2[n] = δ[n-4] + 4δ[n-5] +3δ[n-6] + 2δ[n-7] Write a Matlab program that adds these two signals. Plot the original signals as well as the final result.

-----------------------------TASK 7---------------------------Take a discrete-time signal from . Count the number of samples with amplitude greater than a threshold of 3 and less than a threshold of -3 (use for loop).

-----------------------------TASK 8---------------------------Write your own function to downsample a signal i.e. retain odd numbered samples of the original signal and discard the even-numbered (downsampling by 2). The function must take a signal as input and return the downsampeled version of that signal. See Fig for example. Call this function from a matlab file. your result by using the command “downsample”. Plot the original signal, downsampled signal determined by your program, and downsampled signal obtained by the command downsample.

Fig. DownSampling

-----------------------------TASK 9---------------------------Write your own function to upsample a signal i.e. copy the 1st sample of original signal in the new signal and then place an extra sample of 0, copy the 2nd sample of original signal and then place a 0, and so on. See Fig for example.

Call this function from a matlab file. your result by using the command “upsample”. Plot the original signal, upsampled signal determined by your program, and upsampled signal obtained by the command upsample.

Fig. Upsampling

Lab # 5

OBJECTIVES OF THE LAB

---------------------------------------------------------------------In this lab, we will cover the following topics: • • •

Gain familiarity with Complex Numbers Interpret Phasors & their addition Matlab demo of Phasors

----------------------------------------------------------------------

5.1 COMPLEX NUMBERS A complex number z is an ordered pair (x, y) of real numbers. Complex numbers can be represented in rectangular form (also known as canonical form) as z = x + iy, which is the vector in two-dimensional plane. The horizontal coordinate x is called the real part of z and can be represented as x = Re {z}, while the vertical coordinate y is called the imaginary part of z and represented as y = Imag {z}. That is: z = (x, y) = x + iy = Re {x} + i Imag {x} Another way to represent a complex number is in polar form. In polar form, the vector is defined by its length (r) or magnitude (|z|) and its direction (θ). A rectangular form can be converted into polar form using formulas: |z| = r = (x2 + y2)½ θ = arctan (y/x) z = r ejθ where ejθ = cos θ + i sin θ, and known as the Euler’s formula.

5.2 BUILT-IN MATRIX FUNCTIONS Function Description =============================================== real returns the real part x of z imag returns the imaginary part y of z abs returns the length r of z angle returns the direction θ of z conj returns the complex conjugate ž of z zprint plot vectors in complex z-plane zcat plot vectors in z-plane end-to-end ulot plot a circle with specified center (complex number) and radius ===============================================

Here are some examples: Example To define the complex number, for instance, z = (3, 4) in matlab write in matlab editor >> z = 3 + 4i z= 3.0000 + 4.0000i

Example To find the real and imaginary parts of the complex number, write >> x = real(z) x= 3 >> y = imag(z) y= 4

Example To find the length and direction of z, write >> r = abs(z) r= 5 >> θ = angle(z) θ= 0.9273

Example To find the length and direction of z, write >> zx = conj(z) zx = 3.0000 – 4.0000i

Example To find all the information about a complex number, use the zprint function, i.e. >> zprint(z) Z=

X 3

+

jY 4

Magnitude 5

Phase 0.927

Ph/pi Ph(deg) 0.295

53.13

Example To plot the vector in z-plane, use the zcat function, i.e. >> zcat(z)

Another way to plot is to use the zvect function, which gives the same result as above, i.e. >> z1 = 2 + 3i; >> h = zvect(z1);

Example To plot the circular representation of complex number, ulot function can be used. It takes the radius of circle as first argument, complex number as second argument and any plotting option as third argument. For instance, to draw z1 = 2 + 3i as a dottedgreen circle with radius r = 2 in matlab, write in matlab >> huc = ulot( 1, z1 , ':g');

---------------------------TASK 1-------------------------

Define z1 = -1+j0.3 and z2 = 0.8+j0.7. Enter these in Matlab and plot them with zvect, and print them with zprint.

---------------------------TASK 2-------------------------

Compute the conjugate ź and the inverse 1/z for both z1 and z2 and plot the results. Display the results numerically with zprint.

---------------------------TASK 3-------------------------

Compute z1 +z2 and plot. Use zcat to show the sum as vectors head-to-tail. Use zprint to display the results numerically.

---------------------------TASK 4-------------------------

Compute z1z2 and z1=z2 and plot. Use the zvect plot function to show how the angles of z1 and z2 determine the angles of the product and quotient. Use zprint to display the results numerically.

5.3 COMPLEX EXPONENTIAL SIGNALS The complex exponential signal is defined as x’(t) = A ej(w0t + ø) which is a complex-valued function of t, where the magnitude of x’(t) is |x’(t)| = A

Æ

arg x’(t) = (w0t + ø)

magnitude or length of x’(t) Æ

angle or direction of x’(t)

Using Euler’s formula, it can be expressed in rectangular or Cartesian form, i.e. x’(t) = A ej(w0t + ø) = A cos (w0t + ø) + j A sin (w0t + ø) where A = amplitude, ø=phase shift w0 = frequency in rad/sec

Example clc clear all close all n = 0:1/10:10; k = 5; a = pi/2; x = k * exp(a*n*i); % plot the real part subplot(2,1,1) stem(n, real(x), 'filled') title('Real part of complex exp') xlabel('sample #') ylabel('signal amplitude') grid % plot the imaginary part subplot(2,1,2) stem(n, imag(x), 'filled') title('Imaginary part of complex exp') xlabel('sample #') ylabel('signal amplitude') grid

---------------------------TASK 5------------------------Determine the complex conjugate of the above exponential signal and plot the real and imaginary portions.

---------------------------TASK 6------------------------Generate the complex valued signal y(n) = exp

(-0.1 + j0.3)n

, -10≤n≤10

Plot its magnitude, phase, the real part, and the imaginary part in separate subplots.

---------------------------TASK 7------------------------a) Generate a real-exponential x=an for a=0.7 and n ranging from 0-10. Find the discrete time as well as the continuous time version of this signal. Plot the two signals on same graph (holding both the graphs). b) Repeat the same program with value of a=1.3.

---------------------------TASK 8------------------------Multiply the two discrete signals x1=5exp(i*n*pi/4) and x2= an (use point-by-point multiplication of the two signals). Plot the real as well as the exponential parts for 0 <1 and a rel="nofollow">1.

---------------------------TASK 9-------------------------

Plot the discrete signal x=a^|n| for n ranging from -10 to 10. Draw two subplots for 0 <1 and a rel="nofollow">1.

---------------------------TASK 10-----------------------a) Generate the signal x(t) = Aej(ωt + π) for A = 3, π= -0.4, and ω = 2π(1250). Take a range for t that will cover 2 or 3 periods. b) Plot the real part versus t and the imaginary part versus t. Use subplot(2,1,i) to put both plots in the same window. c) that the real and imaginary parts are sinusoids and that they have the correct frequency, phase, and amplitude.

5.4 ROTATING PHASOR INTERPRETATION The complex exponential signal mentioned above can also be written as x’(t) = A ej(w0t + ø) = A ejø ejw0t = A ejθ(t) ; θ(t) = (w0t + ø)

Æ 1) Æ 2)

Where let X = A ejø in equation 1) denotes the complex amplitude or Phasor. X is the static Phasor as it represents the amplitude and phase shift of the complex exponential signal. While x’(t) in equation 2) denotes the complex number as a vector in complex plane, where the tip of the vector lies on the perimeter of the circle of radius A (amplitude), and it rotates at constant rate with speed w0 (radian frequency). Therefore, x’(t) is known as the rotating Phasor.

Example Matlab Demo of Phasors.

Example – phasor_matlab.m To implement the Rotating Phasor in matlab, consider the following code: A = 1; theta1 = 180/4; t = 0:0.1:2*pi; z = A * exp(j*(t-theta1)); figure; plot(imag(z),real(z),'.'); title('Complex Plane'); xlabel('Real Part'); ylabel('Imag. Part'); axis square; hold on; z1 = []; z2 = []; c = 0; z1(1) = 0 + j*sin(0); z2(1) = 0 + j*sin(0); for i = 0:0.1:2 if(c >= 1) z1(2) = cos((i-0.1)*pi - theta1) + j*sin((i-0.1)*pi - theta1);

plot(z1,'LineWidth', 2, 'Color','w'); z2(2) = cos((i-0.1)*pi - theta1); plot(z2,imag(z(11:12)),'LineWidth', 3, 'Color','w');

end z1(2) = cos(i*pi - theta1) + j*sin(i*pi - theta1); plot(z1,'LineWidth', 2, 'Color','g'); z2(2) = cos(i*pi - theta1); plot(z2,imag(z(11:12)),'LineWidth', 3, 'Color','y'); c = c + 1; pause end grid off

---------------------------TASK 11-----------------------Define z1 = -1+j3 and z2 = 0.8+j1. Enter these in Matlab. Compute z3 = z1 +z2 and plot. Also, draw the rotating Phasor implementation of z1, z2, & z3.

---------------------------TASK 12-----------------------Enhance phasor_matlab.m in such a way that it incorporates the following: a) Generate the real part i.e. Cosine in a separate plot. b) Phasor rotation results into oscillation in corresponding Cosine plot.

Lab # 6

OBJECTIVES OF THE LAB

---------------------------------------------------------------------This lab is mainly concerned with • • • • • •

Generating Sinusoids Sampling a Continuous Time Signal Discrete-Time sinusoids Addition of Sinusoids with Variation in Parameters and their Plots Linear Phase Shift Concept When Dealing With Sum of Sinusoids Three Dimensional Plots

----------------------------------------------------------------------

6.1 GENERATING SINUSOIDS Sinusoidal sequences are implemented using sin() & cos().

Example: a continuous-time sinusoid f0 = 3; A = 5; t = -1:0.005:1; y = A*cos(2*pi*f0*t); figure, plot(t, y,'*:'); xlabel('Time, sec'), ylabel('Amplitude'); title('Graph of sinusoid');

Program: Discrete-Time Sinusoid clc, clear all, close all M=10; n=-3:1/M:3; A=2;

%samples/sec

phase=0; f=1; x=A * sin(2*pi*f*n + phase); stem(n,x,'linewidth', 2) title('Discrete-Time Sine Wave: A sin(2*\pi*f*n + \phi)') xlabel('Time Index') ylabel('Signal Amplitude') axis([n(1) n(end) -A A]) grid

6.2 SAMPLING A CONTINUOUS-TIME SIGNAL A continuous time signal can be sampled using a command: stem(x,y); Following example shows the sampled version of the continuous –time cosine signal.

Example: t = 0:0.0005:1; f = 13; xa = cos(2*pi*f*t); subplot(2,1,1) plot(t,xa);grid xlabel('Time, msec'); ylabel('Amplitude'); title('Continuous-time signal x_{a}(t)'); axis([0 1 -1.2 1.2]) subplot(2,1,2); T = 0.1; n = 0:T:1; xs = cos(2*pi*f*n); k = 0:length(n)-1; stem(k,xs); grid xlabel('Time index n'); ylabel('Amplitude'); title('Discrete-time signal x[n]'); axis([0 (length(n)-1) -1.2 1.2])

---------------------------TASK 1-------------------------

What is the frequency in Hz of the sinusoidal signal? What is the sampling period in seconds?

---------------------------TASK 2------------------------Repeat the program by changing the frequency of the sinusoidal signal to 3 Hz and 7 Hz, respectively. Is there any difference between the corresponding equivalent discrete-time signals. If not, why not?

6.3 ALIASING EFFECT A high frequency gives the sample of lower frequency so that the two can’t be distinguished. If f1 > fs, then its alias are at fa = f1+k*fs;

where k is an integer.

Aliasing exclusively deals with frequencies outside sampling frequency.

Example In this program you will generate a continuous-time equivalent y(t) of the discrete-time signal x[n] to investigate the relation between the frequency of the sinusoidal signal x(t) and the sampling period. To generate the reconstructed signal y(t) from x[n], we x[n] through an ideal low filter. clear all; N=5; fo=3; % maximum frequency in the signal fs=10; % Sampling frequency % Analog signal t=0:0.005:N; x_t=sin(2*pi*fo*t); % Digital signal n=0:1:N*fs; x_n=sin(2*pi*fo*n/fs); j=0; %for k=0:1/fs:N % code for sequence of delayed sinc pulses for k=0:1:N*fs j=j+1; h(j,:)=sinc((t-k/fs)*fs); % Each column represents a delayed sinc end y=x_n*h; plot(n/fs,x_n,'o',t,y);grid; xlabel('Time, msec');ylabel('Amplitude'); title('Reconstructed continuous-time signal y_{a}(t)'); figure;

plot(n/fs,x_n,'o',t,x_t) xlabel('Time, msec');ylabel('Amplitude'); title('continuous-time signal y_{a}(t) sampled');

Example t=0:0.001:1; f1=2*cos(2*pi*1*t); f2=2*cos(2*pi*11*t); n=0:0.1:1; y1=2*cos(2*pi*1*n); y2=2*cos(2*pi*11*n); subplot(2,2,1) plot(t,f1); xlabel('Time'); ylabel('Amplitude'); title('Continous time wave of frequency 1 Hz'); grid; subplot(2,2,2) plot(t,f2); xlabel('Time'); ylabel('Amplitude');

title('Continous time wave of frequency 11 Hz'); grid; subplot(2,2,3) stem(y1); xlabel('sample number'); ylabel('Amplitude'); title('Sampling 1 Hz signal at 10 Hz'); grid; subplot(2,2,4) stem(y2); xlabel('sample number'); ylabel('Amplitude'); title('Sampling 11 Hz signal at 10 Hz'); grid;

---------------------------TASK 3----------------------------

Repeat above program by changing the frequency of the sinusoidal signal to 3 Hz and 7 Hz, respectively. Is there any difference between the corresponding equivalent discrete-time signals and the one generated in Question.

---------------------------TASK 4----------------------------

A square wave can be generated in the same way as you have created a sine wave. Just use the function square instead of sin. Generate a square wave for time -2 to 2 with a frequency of 2Hz.

---------------------------TASK 5---------------------------Generate two 3000 hertz sinusoids with different amplitudes and phases. x1(t) = A1 cos(2π(3000)t + α) x2(t) = A2 cos(2π(3000)t + β)

(a) Select the value of the amplitudes as follows: let A1 = 13 and use your age for A2. For the phases, use the last two digits of your telephone number for α(in degrees), and take β= -30o. When doing computations in Matlab, make sure to convert degrees to radians. (b) Make a plot of both signals over a range of t that will exhibit approximately 3 cycles. Make sure the plot starts at a negative time so that it will include t = 0, and make sure that your have at least 20 samples per period of the wave. (c) that the phase of the two signals x1(t) and x2(t) is correct at t = 0, and also that each one has the correct maximum amplitude. (d) Use subplot (3,1,1) and subplot(3,1,2) to make a three- subplot that puts both of these plots on the same window. (e) Create a third sinusoid as the sum: x3(t) = x1(t) + x2(t). In Matlab this amounts to summing the vectors that hold the samples of each sinusoid. Make a plot of x3(t) over the same range of time as used in the previous two plots. Include this as the third in the window by using subplot (3,1,3). (f) Measure the magnitude and phase of x3(t) directly from the plot. Explain how the magnitude and phase were measured by making annotations on each of the plots.

---------------------------TASK 6---------------------------Generate four sinusoids with the following amplitudes and phases:

i. x1(t) = 5 cos(2π(15)t + 0.5π) ii. x2(t) = 5 cos(2 π (15)t – 0.25π) iii. x3(t) = 5 cos(2 π (15)t + 0.4π) iv. x4(t) = 5 cos(2 π (15)t – 0.9π) (a) Make a plot of all four signals over a range of t that will exhibit approximately 3 cycles. Make sure the plot includes negative time so that the phase at t = 0 can be measured. In order to get a smooth plot make sure that your have at least 20 samples per period of the wave. (b) that the phase of all four signals is correct at t = 0, and also that each one has the correct maximum amplitude. Use subplot (3, 2, i) to make a six- subplot that puts all of these plots on the same page. (c) Create the sum sinusoid via: x5(t) = x1(t) + x2(t) + x3(t) + x4(t). Make a plot of x5(t) over the same range of time as used in the last plot. Include this as the lower in the plot by using subplot (3, 1, 3).

6.4 FOLDING Folding uses the property cos (q) = cos (-q). This causes the frequencies from 0.5 fs to fs become a mirror image of frequencies of 0 to 0.5fs. fapparent = fs – fo, where

fs>fo>0.5fs

Example t=0:0.001:1; f1=2*cos(2*pi*4*t); f2=2*cos(2*pi*6*t); n=0:0.1:1; y1=2*cos(2*pi*4*n); y2=2*cos(2*pi*6*n); subplot(2,2,1) plot(t,f1); xlabel('Time');

ylabel('Amplitude'); title('Continous time wave of frequency 4 Hz'); grid; subplot(2,2,2) plot(t,f2); xlabel('Time'); ylabel('Amplitude'); title('Continous time wave of frequency 6 Hz'); grid; subplot(2,2,3) stem(y1); xlabel('sample number'); ylabel('Amplitude'); title('Sampling 4 Hz signal at 10 Hz'); grid; subplot(2,2,4) stem(y2); xlabel('sample number'); ylabel('Amplitude'); title('Sampling 6 Hz signal at 10 Hz'); grid;

6.5 CREATING PHASE SHIFT Phase shift can be created by adding an angle to 2πft for a sinusoid.

Example clear, close all, clc fs=1000; t=-3:1/fs:3; A=2; phase=0; f=1; x=A * sin(2*pi*f*t + phase); plot(t,x, 'linewidth', 2) title('Continuous-Time Sine Wave: A sin(2*\pi*f*t + \phi)') xlabel('Time Index')

ylabel('Signal Amplitude') axis([t(1) t(end) -A A]) grid

---------------------------TASK 7---------------------------Modify the above program to generate a sine wave with phase shift of +pi/2. Then plot a cosine wave of same frequency, amplitude, and phase shift of 0 in another subplot. Compare both the signals and determine the relationship between the two.

---------------------------TASK 8----------------------------

Write a program to generate a continuous-time sine wave of frequency 3 Hz, positive phase shift of pi/2, and amplitude of 5. Also generate a continuous-time cosine wave of frequency 3 Hz, amplitude of 5, and phase shift of 0. Plot the two signals on separate subplots and properly label them. Determine the relationship between the two signals.

6.6 ADDITION OF SINUSOIDS 6.6.1

CASE 1: When Frequency, Phases, and amplitude of the sinusoids are same t=-2:0.01:2; x1=cos(2*pi*0.5*t); x2=cos(2*pi*0.5*t); x3=x1+x2; subplot(3,1,1); plot(t,x1,'linewidth',3); grid; ylabel('Amplitude'); xlabel('Time'); title('COS WAVE , AMPLITUDE = 1, FREQ = 0.5 HZ, Phase = 0 RADIAN'); subplot(3,1,2); plot(t,x2,'linewidth',3); grid; ylabel('Amplitude'); xlabel('Time'); title('COS WAVE , AMPLITUDE = 1, FREQ = 0.5 HZ, Phase= 0 RADIAN'); subplot(3,1,3); plot(t,x3,'linewidth',3); grid; ylabel('Amplitude'); xlabel('Time'); title('SUM OF THE ABOVE TWO COSINE SIGNALS');

6.6.2 CASE 2: When Frequencies and Phases of the sinusoids are same but Amplitudes are different. t=-2:0.01:2; x1=2*cos(2*pi*0.5*t); x2=cos(2*pi*0.5*t); x3=x1+x2; subplot(3,1,1); plot(t,x1,'linewidth',3); grid; ylabel('Amplitude'); xlabel('Time'); title('COS WAVE , AMPLITUDE = 2, FREQ = 0.5 HZ, Phase = 0 RADIAN'); subplot(3,1,2); plot(t,x2,'linewidth',3); grid; ylabel('Amplitude');

xlabel('Time'); title('COS WAVE , AMPLITUDE = 1, FREQ = 0.5 HZ, Phase= 0 RADIAN'); subplot(3,1,3); plot(t,x3,'linewidth',3); grid; ylabel('Amplitude'); xlabel('Time'); title('SUM OF THE ABOVE TWO COSINE SIGNALS');

6.6.3 CASE 3: When Amplitudes and Phases of the sinusoids are the same but Frequencies are different. t=-2:0.01:2; x1=cos(2*pi*0.5*t); x2=cos(2*pi*1*t); x3=x1+x2; subplot(3,1,1); plot(t,x1,'linewidth',3);

grid; ylabel('Amplitude'); xlabel('Time'); title('COS WAVE , AMPLITUDE = 1, FREQ = 0.5 HZ, Phase = 0 RADIAN'); subplot(3,1,2); plot(t,x2,'linewidth',3); grid; ylabel('Amplitude'); xlabel('Time'); title('COS WAVE , AMPLITUDE = 1, FREQ = 1 HZ, Phase = 0 RADIAN'); subplot(3,1,3); plot(t,x3,'linewidth',3); grid; ylabel('Amplitude'); xlabel('Time'); title('SUM OF THE ABOVE TWO COSINE SIGNALS');

6.6.4 CASE 4: When Amplitudes and Frequencies of the sinusoids are the same but Phases are different t=-2:0.01:2; x1=cos(2*pi*0.5*t); x2=cos((2*pi*0.5*t)+1); x3=x1+x2; subplot(3,1,1); plot(t,x1,'linewidth',3); grid; ylabel('Amplitude'); xlabel('Time'); title('COS WAVE , AMPLITUDE = 1, FREQ = 0.5 HZ, Phase = 0 RADIAN'); subplot(3,1,2); plot(t,x2,'linewidth',3); grid; ylabel('Amplitude'); xlabel('Time'); title('COS WAVE , AMPLITUDE = 1, FREQ = 0.5 HZ, Phase = 1 RADIAN'); subplot(3,1,3); plot(t,x3,'linewidth',3); grid; ylabel('Amplitude'); xlabel('Time'); title('SUM OF THE ABOVE TWO COSINE SIGNALS');

Lab # 7

OBJECTIVES OF THE LAB  

‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐  In this lab, we will cover the following topics:    • Beat Notes  • Amplitude Modulation 

‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ 

7.1 BEAT NOTES When two sinusoidal signals with different frequencies are multiplied, a beat note is produced. To fulfill  the  spectrum  representation,  the  resultant  signal  is  expressed  as  an  additive  linear  combination  of  complex exponential signals. For instance, consider a beat signal as the product of two sinusoids  x(t) = sin(10πt) cos(πt)        = ½ cos(11πt – π/2)  + ½ cos(9πt – π/2)   Or in other words, beat note can be produced by adding two sinusoids with nearly identical frequencies,  i.e.  x(t) = cos(2πf1t)  + cos(2πf2t)  

 

 

 

Æ 1) 

where  f1 = fc – fd 

 

Æ first frequency 

f2 = fc + fd 

 

Æ second frequency 

fc = ½ (f1 + f2)   

Æ center frequency 

fd = ½ (f1 ‐ f2) 

Æ deviation or difference frequency 

 

And, thus   x(t) = sin(2πfdt) cos(2πfct)  

 

 

 

Æ 2) 

 

Æ 3) 

        = ½ cos(2π(fc – fd)t)  + ½ cos(2π(fc + fd)t)      

        = ½ cos(2πf1t)  + ½ cos(2πf2t)  

 

Example To implement the beat note for fc = 200, & fd = 20, consider the following matlab code:   

t = 0:0.00001:0.1;                 fd = 20;                          Ad = 2;                          s1 = Ad * cos(2*pi*fd*t);     

 

% total duration of signal 

     

% deviation or difference frequency  % amplitude of deviation frequency signal  % deviation frequency signal 

fc = 200;                        Ac = 5;                        s2 = Ac * cos(2*pi*fc*t);     

     

% center frequency   % amplitude of center frequency signal  % center frequency signal 

X = s1.*s2;                   

 

% beat signal 

  

    

   % plot deviation frequency signal  figure(1);  plot(t,s1,'linewidth',1.5);  grid;   ylabel('Amplitude');   xlabel('Time');  title('Difference Frequency Signal');    

% plot center frequency signal  figure(2);  plot(t,s2,'linewidth',1.5);   grid;   ylabel('Amplitude');   xlabel('Time');  title('Center Frequency Signal');     % plot the beat signal  figure(3);  plot(t,X,'linewidth',1.5);   grid;   ylabel('Amplitude');   xlabel('Time');  title('Beat Signal'); 

 

Figure: Difference Frequency Signal

 

Figure: Center Frequency Signal

Figure: Beat Signal

--------------------------TASK 1-------------------------Modify the above code in such a way that the resultant beat signal & its envelope both are  shown. 

--------------------------TASK 2-------------------------Write the matlab code that produces a beat signal along with its envelope for frequencies f1 =  191 hertz and f2 = 209 hertz.  

7.2 AMPLITUDE MODULATION   Amplitude modulation is the process of multiplying a low‐frequency signal by a high‐frequency sinusoid.  It is a technique used to broadcast AM radio. The AM signal is the product of the form   x(t) = v(t) cos(2πfct)  

 

 

Æ 4) 

where  the  frequency  of  the  cosine  (fc  hertz)  is  much  higher  than  any  frequencies  contained  in  the  spectrum of v(t), which represent the data signal to be transmitted. The cosine wave is called the carrier  signal, and its frequency is called the carrier frequency.   

Example: Implement the AM signal for the signal given in aforementioned example.  t = 0:0.00001:0.1;            

 

 % total duration of signal 

fd =20;                         Ad = 2;                         s1 = Ad * cos(2*pi*fd*t);    

     

% deviation or difference frequency  % amplitude of deviation frequency signal  % deviation frequency signal 

fc = 200;                       Ac = 5;                         s2 = Ac * cos(2*pi*fc*t);    

     

% center frequency  % amplitude of center frequency signal  % center frequency signal 

B = s1.*s2;                   AM = B + s2;                

   

% beat signal   % AM signal 

  

  

      

  

    

% plot modulating signal  figure(1);  plot(t,s1,'linewidth',1.5);   grid;   ylabel('Amplitude');     xlabel('Time');  title('Modulating Signal');   

 

% plot carrier signal  figure(2);  plot(t,s2,'linewidth',1.5);   grid;   ylabel('Amplitude');   xlabel('Time');  title('Carrier Signal');     % plot the AM signal  figure(3);  plot(t,AM,'linewidth',1.5);   grid;   ylabel('Amplitude');   xlabel('Time');  title('Amplitude Modulated Signal');   

 

Figure: Data Signal

Figure: Carrier Signal

Figure: Amplitude Signal

---------------------------TASK 3------------------------Modify the above code in such a way that the resultant AM signal & its envelope both is shown. 

---------------------------TASK 4------------------------Write the matlab code that produces an AM signal along with its envelope for frequencies f1 =  191 hertz and f2 = 209 hertz.  

---------------------------TASK 5------------------------Write the matlab code that produces an AM signal along with its envelope for the signal                x(t) = ( 5 + 2 cos(40πt) ) cos(400πt). 

---------------------------TASK 6------------------------Write the matlab code that produces an AM signal along with its envelope for the signal                x(t) = ( 12 + 7 sin(πt – π/3) ) cos(13πt).     

              

Lab # 8

OBJECTIVES OF THE LAB  

‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐   This lab aims at the understanding of:   • • • • • • •

Applications of Fourier Series  Synthesis of Square wave  Synthesis of Triangular wave  Synthesis and Analysis of Amplitude Modulated Signals  Chirp Signal  Spectrogram  Frequency Modulated Signal 

‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐       

8.1 FOURIER SERIES   Fourier series theory states that a periodic wave can be represented as a summation of sinusoidal waves  with different frequencies, amplitudes and phase values.  

8.1.1

Synthesis of Square wave

  The square wave for one cycle can be represented mathematically as:  x(t) =           

 

1 

 

0 <= t < T/2 

 ‐1 

 

T/2 <= t < T  

The Complex Amplitude is given by:  Xk   =             

 

 (4/j*pi*k) 

for   

k=±1, ±3, ±5….. 

0 

for   

k=0,±2, ±4, ±6….. 

 

For  f = 1/T = 25Hz, only the frequencies ±25, ±50, ±75 etc. are in the spectrum.   

i.

Effect of Adding Fundamental, third, fifth, and seventh Harmonics

  Example  clear  clc  t=0:0.0001:8;  ff=0.5;  % WE ARE USING SINE FUNCTION BECAUSE FROM EXPONENTIAL FORM OF FOURIER  % SERIES FINALLY WE ARE LEFT WITH SINE   y = (4/pi)*sin(2*pi*ff*t);  % COMPLEX AMPLITUDE = (4/(j*pi*k))  for k = 3:2:7      fh=k*ff;      x = (4/(k*pi))*sin(2*pi*fh*t);      y=y+x;  end 

plot(t,y,'linewidth',1.5);  title('A square wave with harmonics 1st, 3rd, 5th, and 7th');  xlabel('Time');  ylabel('Amplitude');   

   

ii.

Effect of Adding 1st to 17th harmonics

  Example  clear  clc  t=0:0.0001:8;  ff=0.5;  % WE ARE USING SINE FUNCTION BECAUSE FROM EXPONENTIAL FORM OF FOURIER  % SERIES FINALLY WE ARE LEFT WITH SINE   y = (4/pi)*sin(2*pi*ff*t); 

% COMPLEX AMPLITUDE = (4/(j*pi*k))  for k = 3:2:17      fh=k*ff;      x = (4/(k*pi))*sin(2*pi*fh*t);      y=y+x;  end  plot(t,y,'linewidth',1.5);  title('A square wave with harmonics 1st‐17th);  xlabel('Time');  ylabel('Amplitude'); 

   

iii.

Effect of Adding 1st to 27th harmonics

  Example  t=0:0.0001:8; 

ff=0.5;  % WE ARE USING SINE FUNCTION BECAUSE FROM EXPONENTIAL FORM OF FOURIER  % SERIES FINALLY WE ARE LEFT WITH SINE   y = (4/pi)*sin(2*pi*ff*t);  % COMPLEX AMPLITUDE = (4/(j*pi*k))  for k = 3:2:55      fh=k*ff;      x = (4/(k*pi))*sin(2*pi*fh*t);      y=y+x;  end  plot(t,y,'linewidth',1.5);  title('A square wave with harmonics 1st to 27th');  xlabel('Time');  ylabel('Amplitude');   

 

---------------------------TASK 1---------------------------Write a program that plots the signal s(t).   N

sin( 2π nt )

n =1

n

s (t ) = ∑

s (t ) = sin(2π * t ) +

    where n = 1, 3, 5, 7, 9         and N = 9         or 

sin(6π * t ) sin(10π * t ) sin(14π * t ) sin(18π * t ) + + +   3 5 7 9

 

---------------------------TASK 2---------------------------Write a program that plots the signal s(t) in step 2 but with N = 100.     

---------------------------TASK 3---------------------------What do you conclude from TASKS 1 & 2? 

8.2 SYNTHESIS OF A TRIANGULAR WAVE   The Complex Amplitude is given by:  Xk    =    (‐8/*pi^2*k^2)   for k is an odd integer            =    0 

 

 for k for k is an even integer 

For f = 1/T = 25Hz 

  Example: Triangular wave with N=3  t=0:0.001:5;  x=(‐8/(pi*pi))*exp(i*(2*pi*0.5*t));  y=(‐8/(9*pi*pi))*exp(i*(2*pi*0.5*3*t));  s=x+y;  plot(t,real(s),'linewidth',3);  title('Triangular Wave with N=3');  ylabel('Amplitude');  xlabel('Time'); grid; 

 

  Example: Triangular wave with N=11  t=0:0.01:0.25;  ff=25;  x1=(‐8/(pi^2))*exp(i*(2*pi*ff*t));  for k=3:2:21,      fh=ff*k;      x=(‐8/(pi^2*k^2))*exp(i*(2*pi*fh*t));      y=x1+x;  end  plot(t,real(y),'linewidth',3);  title('Triangular Wave with N=11');  ylabel('Amplitude');  xlabel('Time');  grid; 

   

8.3 CHIRP OR LINEARLY SWEPT FREQUENCY   A chirp signal is a sinusoid whose frequency changes linearly from some low value to a high one. To  define the formula for such a signal, following steps can be taken. Since we know that a complex  exponential signal is defined as   x’(t) = A ej (w0t + ø)  And its real part is   x(t)  = Re { A ej (w0t + ø)  }   

 

        = A cos (w0t + ø) 

Then the phase of this signal is the exponent (w0t + ø) that changes linearly with time. The time‐ derivative of phase is w0, which equals the constant frequency.   Thus, the general notation is:  x(t)  = Re { A ej Ψ(t)  }   

 

Where       

        = A cos (Ψ(t))  Ψ(t) Æ represents the phase as a function of time  Ψ(t) = 2πµt2 + 2πf0t + ø  

  The derivative of Ψ(t) yields an instantaneous frequency that changes linearly versus time.      fi (t) = 2µt + f0     The slope of fi(t) is equal to 2Ψ and its intercept is equal to f0. If the signal starts at t = 0, then f0 is also  the starting frequency. The frequency variation produced by the time‐varying phase is called frequency  modulation, and this class of signals is called FM signals. Finally, since the linear variation of the  frequency can produce an audible sound similar to a siren or a chirp, the linear‐FM signals are also called  chirps." 

Example: The following Matlab code synthesizes a chirp: fsamp = 8000;        dt = 1/fsamp;        dur = 1.8;        tt = 0 : dt : dur;        psi = 100 + 2*pi*(200*tt + 500*tt.*tt);  A = 7.7;         xx = real( A * exp (j*psi) );    sound( xx, fsamp );       

% sampling frequency  % increment value  % total duration    % time vector    % instantaneous phase: ø = 100, f0 = 200, µ = 500   % Amplitude  % chirp signal  % play the signal at given sampling frequency   

---------------------------TASK 4---------------------------(a) Determine the range of frequencies (in hertz) that will be synthesized by above mentioned  Matlab script. Make a sketch by hand of the instantaneous frequency versus time. What are the  minimum and maximum frequencies that will be heard" signal according to the following comments:  function xx = mychirp( f1, f2, dur, fsamp )    % MYCHIRP generate a linear‐FM chirp signal  % usage: xx = mychirp( f1, f2, dur, fsamp )  % f1 = starting frequency  % f2 = ending frequency  % dur = total time duration  % fsamp = sampling frequency (OPTIONAL: default is 8000)    if( nargin < 4 ) %‐‐ Allow optional input argument  fsamp = 8000;  end 

 

(c) Generate a chirp sound to match the frequency range of the chirp in part (a). Listen to  the chirp using the sound function. Also, compute the spectrogram of your chirp using the  Matlab function: specgram(xx,[],fsamp). 

---------------------------TASK 5---------------------------Use your Matlab function mychirp to synthesize a “chirp" signal for your lab report. Use the  following parameters:    1. A total time duration of 3 secs with a D/A conversion rate of fs = 8000 Hz.  2. The instantaneous frequency starts at 15,000 Hz and ends at 300 Hz.    Listen to the signal. What comments can you make regarding the sound of the chirp (e.g. is it  linear)? Does it chirp down, or chirp up, or both? Create a spectrogram of your chirp signal. Use  the sampling theorem (from Chapter 4 in the text) to help explain what you hear and see.     

8.4 FREQUENCY MODULATED SIGNAL   In the constant frequency sinusoid, the argument of the cosine is also the exponent of the  complex exponential, so the phase of this signal is the exponent (2Ψf0t + Ψ). This phase  function changes linearly versus time, and its time derivative is 2Ψf0 which equals the constant  frequency of the cosine.    A generalization is available if we adopt the following notation for the class of signals with time‐ varying phase  (3)

The time derivative of the phase from (3) gives a frequency 

but we prefer units of hertz, so we divide by 2Ψ to define the instantaneous frequency: 

8.5 SPECTROGRAM   It is often useful to think of signals in  of their spectra. A signal's spectrum is a  representation of the frequencies present in the signal. For a constant frequency sinusoid the  spectrum consists of two spikes, one at 2Ψf0, the other at ¡2Ψf0. For more complicated signals  the spectra may be very interesting and, in the case of FM, the spectrum is considered to be 

time‐varying. One way to represent the time‐varying spectrum of a signal is the spectrogram  (see Chapter 3 in the text). A spectrogram is found by estimating the frequency content in short  sections of the signal. The magnitude of the spectrum over individual sections is plotted as  intensity or color on a two‐dimensional plot versus frequency and time.    There are a few important things to know about spectrograms:    1. In Matlab, the function spectrogram computes the spectrogram. Type help spectrogram  to learn more about this function and its arguments.    2. Spectrograms are numerically calculated and only provide an estimate of the time‐ varying frequency content of a signal. There are theoretical limits on how well they can  actually represent the frequency content of a signal. This problem will be treated in the  future lab when we use the spectrogram to extract the frequencies of piano notes.    Beat notes provide an interesting way to investigate the time‐frequency characteristics of  spectrograms.    Although some of the mathematical details are beyond the reach of this course, it is not  difficult to understand the following issue: there is a fundamental trade‐o® between knowing  which frequencies are present in a signal (or its spectrum) and knowing how those frequencies  vary with time. As mentioned previously in Section 1.4, a spectrogram estimates the frequency  content over short sections of the signal. Long sections give excellent frequency resolution, but  fail to track frequency changes well. Shorter sections have poor frequency resolution, but good  tracking. This trade‐off between the section length (in time) and frequency resolution is  equivalent to Heisenburg's Uncertainty Principle in physics.     A beat note signal may be viewed as a single frequency signal whose amplitude varies with  time, or as two signals with different constant frequencies. Both views will be useful in  evaluating the effect of window length when finding the spectrogram of a beat signal.   

---------------------------TASK 6--------------------------- (a) Create and plot a beat signal with 

(i) f¢ = 32 Hz  (ii) Tdur = 0:26 sec  (iii) fs = 8000 Hz, or 11,025 Hz  (iv) f0 = 2000 Hz   

(b) Find the spectrogram using a window length of 2048 using the commands:  specgram(x,2048,fsamp); colormap(1‐gray(256)).  Comment on what you see.    (c) Find the spectrogram using a window length of 16 using the commands:  specgram(x,16,fsamp); colormap(1‐gray(256)). Comment on what you see. 

   

Lab # 9

OBJECTIVES OF THE LAB  

‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐  This lab aims at the understanding of:   • • • • •

Generating unit impulse and unit step sequences  Basic signal operations  Characterization of LSI systems  Implementation of Running Average Filter (Causal and Non‐Causal)  Delay Filter 

‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐       

9.1     GENERATING UNIT IMPULSE AND UNIT STEP SEQUENCES   Use matlab commands zeros and ones. 

Example: Unit Impulse Sequence    n=‐10:10;    % unit impulse  x1=[zeros(1,10) 1 zeros(1,10)];    stem(n,x1,'filled');  xlabel('sample #');  ylabel('signal amplitude');  title('Unit impulse');  axis([‐10 10 ‐1 2]); 

 

Example: Unit Step Sequence    n= ‐10:10;    %unit step  x1=[zeros(1,10) ones(1,11)];    stem(n,x1,'filled');  xlabel('sample #');  ylabel('signal amplitude');  title('Unit step');  axis([‐10 10 ‐1 2]); 

   

---------------------------TASK 1------------------------Plot the signum sequence. It can be defined as:  Sign(n)=

1, -1, 0,

for n>0 for n<0 for n=0

9.2     BASIC SIGNAL OPERATIONS: 1) Signal Shifting clc  clear all  close all    n=0:0.002:4;  x=sin(2*pi*1*n);    subplot(2,1,1);  plot(n,x,'linewidth',2);  title('Original Signal');  xlabel('Time');  ylabel('Signal Amplitude');  axis([‐3 3 ‐1 1]);  grid;    subplot(2,1,2);  plot(n‐1,x,'linewidth',2);  title('Advanced signal');  xlabel('Time');  ylabel('Signal Amplitude');  axis([‐3 3 ‐1 1]);  grid; 

 

---------------------------TASK 2------------------------Delay the above signal by 1 sec. Plot both the delayed & original signal on the same figure. 

2) Signal Flipping clear  n=‐1:1/1000:1;  x1=5*sin(2*pi*1*n);    subplot(2,1,1);  plot(n,x1, 'g', 'linewidth',2);  axis([‐1 1 ‐5 5]);  xlabel('time');  ylabel('signal amplitude');  title('Original sine wave');  grid;   

subplot(2,1,2);  plot(‐n,x1, 'r', 'linewidth',2);  axis([‐1 1 ‐5 5]);  xlabel('time');  ylabel('signal amplitude');  title('Flipped Sine Wave');  grid; 

 

---------------------------TASK 3------------------------Flip the following signal: 

Y= 5exp (i*n*pi/4) Plot the original signal as well as the flipped one in the same figure. 

---------------------------TASK 4------------------------Flip the following signal:  X[n]= 2δ[n] +  5δ[n‐1] + 8δ[n‐2] + 4δ[n‐3] + 3δ[n‐4]  Plot the original signal as well as the flipped one in the same figure.  3) Amplitude Scaling clear  n=1:7;  x=[1 2 2 3 2 2 1];    subplot(2,1,1);  stem(n,x, 'filled');  title('Original signal');  xlabel('Time index');  ylabel('Signal Amplitude');  axis([1 7 0 4]);  grid;    S=2;  subplot(2,1,2);  stem(n,S*x, 'filled');  title('Amplitude Scaled signal');  xlabel('Time index');  ylabel('Signal Amplitude');  axis([1 7 0 8]);  grid;     

 

---------------------------TASK 5------------------------Scale the continuous-time sinusoid of program 3 (i.e. the one used in signal shifting operation) by a factor of 2. 4) Time Scaling %Decimation (down‐sampling)  clear    n=‐2:1/1000:2;  x1=sin(2*pi*2*n);  x2=decimate(x1,2);    subplot(2,1,1);  plot(x1);  title('Original signal');  xlabel('Sample Number'); 

ylabel('Signal Amplitude');  axis([0 4000 ‐1 1]);  grid;    subplot(2,1,2);  plot(x2);  title('Decimated signal');  xlabel('Sample Number');  ylabel('Signal Amplitude');  axis([0 2000 ‐1 1]);  grid; 

   

---------------------------TASK 6------------------------Use interp command in the above program to interpolate (up-sample) the signal by a factor of 2.

5) Amplitude Clipping clear  x=[3 4 4 2 1 ‐4 4 ‐2];  len=length(x);  y=x;  hi=3;  lo=‐3;  for i=1:len  if(y(i)>hi)       

y(i)=hi;  elseif(y(i)
end    subplot(2,1,1);  stem(x,'filled');  title('original signal');  xlabel('Sample number');  ylabel('Signal Amplitude');    subplot(2,1,2);  stem(y,'filled');  title('Clipped Signal');  xlabel('Sample number');  ylabel('Signal Amplitude');

 

6) Signal Replication clear  x=[1 2 3 2 1];  y=[x x x x];    subplot(2,1,1);  stem(x,'filled');  title('Original Signal');  xlabel('Sample Number');  ylabel('Signal Amplitude');  axis([1 20 0 3]);  grid;    subplot(2,1,2);  stem(y,'filled');  title('Replicated Signal'); 

xlabel('Sample Number');  ylabel('Signal Amplitude');  axis([1 20 0 3]);  grid; 

 

9.3 CHARACTERIZATION OF LSI SYSTEMS   There are different ways to characterize discrete‐time LSI systems including difference equation,  impulse response, transfer function, and frequency response. We are considering here only difference  equation. 

Difference Equation A discrete‐time LSI system is characterized by the following general form difference equation  y(n) + a1y1(n-1) + a2y2(n-2) + … + aNyN(n-N) = b0x0(n) + b1x1(n-1) + …+ bMxM(n-M) where     ¾ x(n) and y(n) are input and output, respectively.   ¾ Finite integers M and N represent the maximum delays in the input and output respectively. 

¾ The constants ai and bi are called the filter coefficients.  Difference equation can be implemented in Matlab as follow:  y = filter (b,a,x) which takes the input signal in vector x and filter coefficient a and b. The Vector B constitutes the feed‐ forward filter coefficients and vector A constitutes the feed‐backward filter coefficients. In case of FIR  filters, the feed‐backward filter coefficient consists of a single 1. 

Example  y (n) + 2y(n‐1) + 3y(n‐2) = x(n) + 3x(n‐1) + x(n‐2)  A = coefficients of y(n) = [1 2 3]  B = coefficients of x(n) = [1 3 1]  X = input sequence = [1 2 3 4 5 6 ] 

Program  A = [1 2 3];   B = [1 3 1];   X = [1 2 3 4 5 6];   Y =filter (B, A, X) 

Output :   

 

Y = 

 

 

     1     3    1     4     9    ‐5 

---------------------------TASK 7------------------------Given the following difference equation and input signal, calculate the output.  ¾ 5y(n) = ‐4x(n) + x(n‐1) + 6x(n‐3)  ¾ X = input sequence = [1 3 5 7 9]     

9.4 FILTERING SIGNALS   Matlab provides the command Filter for developing One‐dimensional digital filter.   Y = filter(B,A,X)   filters the data in vector X with the filter described by vectors A and B to create the filtered data Y.  The  filter is a "Direct Form II Transposed" implementation of the standard difference equation:     a(1)*y(n) =   b(1)*x(n) + b(2)*x(n‐1) + ... + b(nb+1)*x(n‐nb)  ‐ a(2)*y(n‐1) ‐ ... ‐ a(na+1)*y(n‐na) 

For FIR filters, keep a=1;   

Example    clc  clear    b = [1 2 3 4 5 4 3 2 1];  a = 1;  x = cos(0.2*pi*[0:20]);  y = filter(b,a, x);    figure;  subplot(2,1,1);  stem (x);  title('Discrete Filter Input x[n]');  xlabel('index, n');  ylabel('Value, x[n]');  axis([0 21 ‐1 1]);  grid;    subplot(2,1,2);  stem(y, 'r');  title('Discrete Filter Output y[n]');  xlabel('index, n');  ylabel('Value, y[n]');  axis([0 21 ‐20 20]);  grid; 

   

9.4.1

3-Point Averaging filter clear  n=‐3:7;    % input signal  x=[0 0 0 2 4 6 4 2 0 0 0];    % filter coefficients  b=[1/3  1/3  1/3];       

 

 

%feedforward filter coefficeints 

a=1;                      

 

 

 % filter coefficients 

  % output signal  y=filter(b,a,x);     

figure;  subplot(2,1,1);  stem(n,x, 'filled');  xlabel('sample number');  ylabel('signal amplitude');  title('origianl signal');  grid;    subplot(2,1,2);  stem(n,y, 'filled');  xlabel('sample number');  ylabel('signal amplitude');  title('Result of 3‐point averager: y[n]=(x[n]+x[n‐1]+x[n‐2])/3');  grid; 

   

---------------------------TASK 8------------------------Modify above program to implement:    1. A 3‐point non‐causal averager  2. A centralized average 

---------------------------TASK 9------------------------A filter has coefficients:  bk={1  3  ‐2   4  1}  Input signal  x[n]= δ[n] ‐  3δ[n‐1] + 3δ[n‐2] + 4δ[n‐3] ‐ 2δ[n‐4]    is applied. Find the output y[n]. 

---------------------------TASK 10-----------------------A non‐recursive discrete‐time system has filter coefficients {1, ‐2, 3, 2}. An input signal (shown  in Fig) is applied to the system. Determine the discrete‐time output of the filtering. Write code  to plot the original signal as well as the output on two subplots of the same figure. 

 

2 2 1

 

1

1

  Example 

‐2  ‐1    0   1   2   3     4     5     6    7

n=0:40;  % input signal  x=(1.02).^n + 0.5*cos(2*pi*n/8 + pi/4);    % filter coefficients  b1=[1/3 1/3 1/3];  a1=1;    % plot of original signal  figure; 

subplot(2,1,1);  stem(n,x, 'filled');  title('input signal');  xlabel('sample number');  ylabel('signal amplitude');  grid;    % output of 3‐ point filtering  y1=filter(b1,a1,x);    subplot(2,1,2);  stem(n,y1, 'filled');  title('output of 3‐point filtering');  xlabel('sample number');  ylabel('signal amplitude');  grid; 

---------------------------TASK 11-----------------------Develop a causal FIR filter that averages five samples of the input signal  x[n]= 2δ[n] +  5δ[n‐1] + 8δ[n‐2] + 4δ[n‐3] + 3δ[n‐4]   

---------------------------TASK 12-----------------------Apply 7‐point filtering on the following signal:   x=(1.02).^n + 0.5*cos(2*pi*n/8 + pi/4);   

0≤n≤40 

9.5 DELAY FILTER A delay system is a simple FIR filter. A delay by 2 is given by filter coefficients {0, 0, 1}. Its  implementation on input signal i.e. x[n] = {2 4 6 4 2} is shown in figure below: 

 

---------------------------TASK 13-----------------------Given the signal   x[n]= δ[n] +  4δ[n‐1] +3δ[n‐2] + 2δ[n‐3]  Delay this signal by 2 using the above mentioned filter coefficients by means of filter command.  

---------------------------TASK 14-----------------------A delay filter delays an input signal. The amount of delay depends on filter coefficients. Design a  delay filter that delays an input signal by 6 units.  

---------------------------TASK 15-----------------------Use the above filter to delay the input signal:  x[n]= δ[n] - 4δ[n-1] +3δ[n-2] + 2δ[n-3] - 6δ[n-4] + 2δ[n-6] Write code to plot the input signal as well as the output signal.  

---------------------------TASK 16-----------------------Create a function that takes a signal and a delay count as arguments. The purpose of the  function is to delay the signal by that amount (use filter command). Call this function from a  main program to delay a required signal. 

Lab # 10

OBJECTIVES OF THE LAB  

‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐  This lab aims at the understanding of:   • • • • • • •

Making Signals Causal and Non‐Causal  Convolution  Properties of Convolution  Implementation of Low  FIR Filters and its Frequency Response  Implementation of High  FIR Filters and its Frequency Response  Implementation and Frequency Response of a Moving Average Filter  Implementation of Band  Filters 

 

‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐   

10.1  

MAKING SIGNALS CAUSAL AND NON-CAUSAL

Causal  Signal:  A  signal  is  said  to  be  causal  if  it  is  zero  for  time<0.  A  signal  can  be  made  causal  by  multiplying it with unit step.  

  Example  clc  clear all  close all     t = ‐2:1/1000:2;  x1 = sin(2*pi*2*t);     subplot(3,1,1);  plot(t,x1,'LineWidth',2);  xlabel('time');  ylabel('signal amplitude');  title('sin(2*\pi*f*t)');     u = (t>=0);  x2 = x1.*u;     subplot(3,1,2);  plot(t,u, 'r','LineWidth',2);  xlabel('time');  ylabel('Signal Amplitude');  title('Unit Step');     subplot(3,1,3);  plot(t,x2, 'k','LineWidth',2);  xlabel('time');  ylabel('signal amplitude');  title('causal version of sin(2*\pi*f*t)');     figure;  plot(t,x1,t,u,'‐.',t,x2,'LineWidth',2);  text(0,1.2,'u(t)','FontSize',16);  text(‐1.2,‐1.1,'x(t)','FontSize',16);  text(0.8,‐1.1,'x(t)*u(t)','FontSize',16);  axis([‐2 2 ‐1.5 1.5]); 

   

 

--------------------------TASK 01------------------------Sample the above signal to get its discrete‐time counterpart (take 10 samples/sec). Make the  resultant signal causal. Display the lollypop plot of each signal. 

--------------------------TASK 02------------------------A signal is said to be anti‐causal if it exists for values of n<0. Make the signal given in above  example anti‐causal.   

--------------------------TASK 03------------------------Create a function by name of sig_causal in matlab that has two input arguments: (i) a discrete‐ time signal, and (ii) a position vector. The function should make the given signal causal and  return the resultant signal to the calling program.   A non‐causal signal is shown in the Fig below. Write matlab code to make the signal causal by  calling the above‐mentioned function. Plot the original non‐causal signal and the resultant  causal signal.   

2 2

 

1

 

1

1

 

 

‐2  ‐1    0   1   2   3     4     5     6    7  

10.2 CONVOLUTION: Use the matlab command Conv(h,x) to find convolution where  H = impulse response   X = input signal 

Example  clc  clear all  close all 

  h = [1 2 3 4 5 4 3 2 1];  x = sin(0.2*pi*[0:20]);  y = conv(h, x);    figure(1);  stem (x);  title('Discrete Filter Input x[n]');  xlabel('index, n');  ylabel('Value, x[n]');    figure (2);  stem(y, 'r');  title('Discrete Filter Output y[n]');  xlabel('index, n');  ylabel('Value, y[n]');

 

  Even though there are only 21 points in the x array, the conv function produces 8 more points because it  uses the convolution summation and assumes that x[n] = 0 when n>20.   

--------------------------TASK 04------------------------Convolve the following signals:  x = [2 4 6 4 2];  h = [3 ‐1 2 1];  plot the input signal as well as the output signal.   

--------------------------TASK 05------------------------Convolve the signal x[n]=[1 2 3 4 5 6] with an impulse delayed by two samples. Plot the original  signal and the result of convolution.  

 

--------------------------TASK 06------------------------Convolution is associative.  Given the three signal x1[n], x2[n], and x3[n] as:  x1[n]= [ 3 1 1]  x2[n]= [ 4 2 1]  x3[n]=[ 3 2 1 2 3]  Show that (x1[n] * x2[n]) * x3[n] = x1[n] * (x2[n] * x3[n])   

--------------------------TASK 07------------------------Convolution is commutative. Given x[n] and h[n] as:  X[n]=[1 3 2 1]  H[n]=[1 1 2]  Show that x[n] * h[n] = h[n] * x[n]   

--------------------------TASK 08------------------------Determine h(n) for the system  10

y (n) = ∑ kx(n − k )   k =0

When x(n)= 2δ[n]. Plot the input, impulse, and the output. 

 

 

 

--------------------------TASK 09------------------------Given the impulse response of a system as:  h[n] =  2δ[n] +  δ[n‐1] + 2δ[n‐2] + 4δ[n‐3] + 3δ[n‐4]  If the input  x[n] =  δ[n]  +  4δ[n‐1] +3δ[n‐2] + 2δ[n‐3] is applied to the system, determine the  output of the system. 

 

--------------------------TASK 10------------------------Two systems are connected in cascade.     y[n]      x[n]         h1[n]         h2[n]      h1[n]=[1 3 2 1]  h2[n]=[1 1 2]  If the input x[n] =  δ[n] +  4δ[n‐1] +3δ[n‐2] + 2δ[n‐3] is applied, determine the output.   

--------------------------TASK 11------------------------Given the signals:     

x1[n]= 2δ[n] ‐  3δ[n‐1] + 3δ[n‐2] + 4δ[n‐3] ‐ 2δ[n‐4]  x2[n]= 4δ[n] + 2δ[n‐1] + 3δ[n‐2] ‐ δ[n‐3] ‐ 2δ[n‐4]  x3[n]= 3δ[n] + 5δ[n‐1]‐ 3δ[n‐2] +4 δ[n‐3]  

 that  

 

•

x1[n] * (x2[n] * x3[n]) = (x1[n] * x2[n]) * x3[n] 

•

x1[n] * x2[n]= x2[n] * x1[n] 

Lab # 11

OBJECTIVES OF THE LAB  

‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐  This lab aims at the understanding of:   • • • • •

Implementation of Low  FIR Filters and its Frequency Response  Implementation of High  FIR Filters and its Frequency Response  Implementation and Frequency Response of a Moving Average Filter  Implementation of Band  Filters  Deg filters using matlab command FIR1 

 

‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐   

 

11.1 FREQUENCY RESPONSE 11.1.1 Using Freqz   Given the coefficients of a filter, frequency response can be determined for a particular range of  frequencies by using “Freqz”. The syntax is:  H = freqz(b, a, w)  Where b and a are the filter coefficients and w is the range of frequencies for which the  frequency response is desired to be determined. To plot the magnitude and phase of frequency  response, use abs and angle respectively.  

   

Example: A simple Low filter  clc  clear all  close all    b = [1 2 1];  a = 1;  w = ‐3*pi:1/100:3*pi;    H = freqz(b, a, w);    subplot(3,1,1);  plot(w, abs(H), 'LineWidth', 2);  title('Magnitude of Frequency Response of filter with coefficients b_k={1, 2, 1}');  xlabel('\omega');  ylabel('H(\omega)');  grid;    subplot(3,1,2);  plot(w, abs(H), 'r ', 'LineWidth', 2); 

title('Zoomed view of the above graph from ‐\pi to \pi');  xlabel('\omega');  ylabel('H(\omega)');  axis([‐pi pi min(H) max(H)+0.5]);  grid;    subplot(3,1,3);  plot(w, angle(H), 'k ', 'LineWidth', 2);  title('Phase of Frequency Response of filter with coefficients b_k={1, 2, 1}');  xlabel('\omega');  ylabel('H(\omega)');  grid; 

--------------------------TASK 01------------------------Plot  the  magnitude  and  phase  of  the  frequency  response  of  an  FIR  filter  with  coefficients                  b = {1 ‐2 4 ‐2 1}. Determine the type of filter.  

--------------------------TASK 02------------------------Repeat task 9 for filter coefficients b={1, ‐1}. 

--------------------------TASK 03------------------------A causal three point running average filter is given by   Y[n]=(x[n] + x[n‐1] + x[n‐2])/3;  Find the frequency response of this filter from ‐3π to +3π. 

Note: Freqz command can also be used as:  [H W] = freqz(b, a, n);  Where H contains the frequency response, W contains the frequencies between 0‐π where the  response is calculated, and n is the number of points at which to determine the frequency  response. If n is missing, this value defaults to 512. 

   

Example: A simple Low filter  clc  clear all  close all    b = [1 2 3 2 1];  a = 1;    [H W] = freqz(b, a);    subplot(2,1,1);  plot(W, abs(H), 'linewidth', 1.5);  title('Magnitude of Frequency Response of filter with coefficients b_k={1, 2, 3, 2, 1}');  xlabel('\omega ‐‐‐‐>');  ylabel('H(\omega) ‐‐‐‐>'); 

axis tight;  grid;    subplot(2,1,2);  plot(W, angle(H),'g', 'linewidth', 1.5);  title('Phase of Frequency Response of filter with coefficients b_k={1, 2, 3, 2, 1}');  xlabel('\omega ‐‐‐‐>');  ylabel('H(\omega) ‐‐‐‐>');  axis tight;  grid; 

 

--------------------------TASK 04------------------------Using  the  above‐mentioned  syntax,  determine  the  frequency  response  of  a  filter  with  filter  coefficients {‐1, 3, ‐1}. Take 1024 points. 

11.2 FIR1 Fir1 is used to design filter coefficients. The syntax of this command is:  b = fir1(N, Wc)  Where b contains the filter coefficients returned by this command, N is the order of filter, and wc  is the cutoff frequency normalized between 0 ‐‐ 1 where 1 refers to π. The number of coefficients  created is one more than the order of filter. By default, a low filter is created. A high  filter can be created using the argument ‘high’ at the end in the above command. A band  filter can be created using [Wa Wb] instead of Wc where frequencies between Wa and Wb are  allowed to . If such a vector is specified instead of Wc, a band filter is created by default.  A bandstop filter can be created by using the argument ‘stop’ in the aforementioned command.  B=fir1(N, Wc)  

 

 

//creates coefficients for low filter 

B=fir1(N, Wc, ‘high’) 

 

//creates coefficients for high filter 

B=fir1(N, [Wa Wb]) 

 

//creates coefficients for band filter 

B=fir1(N, [Wa Wb], ‘stop’) 

 //creates coefficients for bandstop filter 

After the coefficients are obtained, freqz command can be used to determine the frequency  response.  Note: For High and bandstop filters, N must be even. Else Matlab will increase the order  itself by 1. 

   

Example: Deg High  filter  % High  filter  clc, clear all, close all    w = ‐pi:1/100:pi;  b = fir1(4, 1/3, 'high')  H = freqz(b, 1, w);  subplot(2,1,1);  plot(w/pi, abs(H) , 'linewidth', 2); 

title('High  filter with cut‐off frequency of w_c=\pi/3');  xlabel('Normalized \omega = x\pi');  ylabel('H(\omega)');  axis tight;  grid;    b = fir1(1000, 1/3, 'high');  H = freqz(b, 1, w);    subplot(2,1,2);  plot(w/pi, abs(H), 'g', 'linewidth', 2);  title('High  filter with cut‐off frequency of w_c=\pi/3');  xlabel('Normalized \omega = x\pi');  ylabel('H(\omega)');  grid; 

--------------------------TASK 05------------------------Design a low‐ filter with 7 coefficients and cutoff frequency of π/2. 

--------------------------TASK 06------------------------Design a band‐ filter which allows only the frequencies between π/6 and 2π/3 to  and  blocks the rest of frequencies. Design it for order 10 and 1000.   

--------------------------TASK 07------------------------Design a band‐stop filter which stops the frequencies between π/3 and π/2 and allows the rest  of frequencies to . Design it for order 10 and 1000. 

 

--------------------------TASK 08------------------------Magnitude of an ideal high  filter can be extracted from an ideal low filter as        Hhp=1‐Hlp and vice‐versa. Design a high‐ filter with order 1000 and Wc=π/4. Using the  high‐ filter’s frequency response, determine the frequency response of a low  filter  with Wc=π/4. (Hint: use the values 1‐abs(H) ).  

 

--------------------------TASK 09------------------------Magnitude of an ideal bandstop filter can be extracted from a band filter as Hbs=1‐Hbp and  vice versa. Create a band filter with order 1000 which allows only the frequencies between  π/6 and 2π/3 to  and blocks the rest of frequencies. Create a bandstop filter from this  band filter. 

--------------------------TASK 10------------------------Impz command is used to find the impulse response of any filter if the coefficients of that filter  are known. The syntax of this command is:  h=impz(b, a);  or [h t]=impz(b,a); where t is the vector of time.   or [h t]= impz(b, a, n); where n is the number of samples of the impulse response. 

Determine the impulse response of an FIR filter with coefficients {2, 3, 1, 3, 2}.  In the second  part, determine 10 samples of the above‐mentioned impulse response.   

--------------------------TASK 11------------------------ (a)  Design a band‐ filter of order 15 that allows only a band of frequencies between π/6  and π/3 and suppresses all the other frequencies. Determine the gain of this system in the  range ‐3π to +3π.  (b) Design this band‐ filter for an approximation to ideal filtering.  

--------------------------TASK 12------------------------Write a program to find the impulse response of a 11‐point running average filter.  

--------------------------TASK 13------------------------Design a high‐ filter with 10 coefficients. The cutoff frequency of this filter should be π/3.  Determine the frequency response of the filter from ‐4π to +4π. Plot the magnitude and phase  response of this filter.   

         

Lab # 12

OBJECTIVES OF THE LAB  

‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐  This lab aims at the understanding of:   • • • •

Z‐Plane and Unit Circle  Pole Zero Plots  Implementation and Frequency Response of Nulling Filter   Implementation and Frequency Response of Running Sum Filter and its Pole Zero Plot 

‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐ 

     

12.1 THE THREE DOMAINS   There are three domains in which a signal can be analyzed: 

n‐Domain:     

It deals with sequences, impulse response and difference equation. 

Ŵ‐Domain:       

Frequency domain deals with spectrum and frequency response. 

Z‐Domain:  • • •

 

This domain deals with the Z‐operator, poles and zeros.  There are the issues of stability and it will eventually lead us to IIR filters.  Frequency transform is also considered as a special case of z‐transform. 

  Definition of Z‐Transform   Consider a finite length signal x[n]: 

   

 

x[n] =k=0∑M x[k]δ[n‐k] 

The Z‐Transform of such a signal is defined by: 

   

   

X(z) = k=0∑M x[k] z‐k 

 

System Function: 

 

 

H(z) = k=0∑M x[k] z‐k  

Where, 

    Y(z) = H(z)X(z)    The Z Plane And Unit Circle  • • •

The z operator is complex exponential of the form z = rejw   As ‘r’ becomes 1 the z‐transform reduces to Fourier transform.  This forms a circle of radius 1 and it is called unit circle.   

Poles:   

These are the values of z on which the value of transform is undefined. 

Zeros:     

 

These are the value of z on which the value of transform is zero. 

Example:  Consider the system function: 

 

 

H(z) = 1‐ 2z‐1 + 2 z‐2 ‐ z‐3 

  Matlab Command: 

 

 

zplane(b,a) 

where ,b and a are row vectors containing transfer function polynomial coefficients  b =  Coefficients of Numerator  a =  Coefficients of Denominator    We can rewrite the system function as: 

           H(z) = (z3 ‐2 z2 +2 z ‐1) / z3                      = (z‐1)(z‐ ej(pi/3)n) (z‐ e‐j(pi/3)n)/ z3    Matlab Command: 

                           zplane(Z,P)    Plots the zeros Z and poles P (in column vectors) with the unit circle for reference  where ,  Z = Zeros of System Function  P = Poles of System Function   

  Example: Pole Zero Plot of a System Function with Coefficients Given By B    b = [1 ‐2 2 ‐1];  a = [1];    figure;  zplane(b,a);  title('Pole‐zero plot of H(z) = 1 ‐ 2z(‐1) + 2z(‐2) ‐z(‐3)');   

% FREQUENCY RESPONSE OF THE SYSTEM  w = ‐pi:pi/500:pi;  h = freqz(b,a,w);    figure;  plot(w,abs(h),'linewidth',2.5);  title('Magnitude of frequency response');  xlabel('Normalized radian Frequency');  ylabel('Amplitude') 

 

   

Example: % Running Sum Filter Quite Similar To Moving Average  b = ones(1,11);% coefficients of numerator of z transform  a = [1];% coefficients of denominator of z transform    figure(1);  zplane(b,a);  title('Pole‐zero plot of 11‐Point Running sum filter');    w = ‐pi:pi/500:pi;  h = freqz(b,a,w);  figure;  plot(w,abs(h),'linewidth',2.5);  title('Magnitude of frequency response 11‐Point Running sum filter');  xlabel('Normalized radian Frequency');  ylabel('Amplitude') 

 

 

--------------------------TASK 01------------------------Design a Band  Filter using the running sum filter.  Note: Take the length of the filter of your own choice. 

12.2 STABILITY OF A SYSTEM   A system is said to be stable if all the poles lie inside the unit circle. This can be determined by pole‐zero  plot of a filter.  

  Example  b = [5 3 1 2 1 3 5] %coefficients of numerator of z transform  a = [1];% coefficients of denominator of z transform    figure(1);  zplane(b,a);  title('Pole‐zero plot of an FIR filter'); 

 

Since the poles lie inside the unit circle, the system is stable. 

--------------------------TASK 02------------------------Determine the frequency response as well as the pole‐zero plot for an FIR filter determined by  the feed‐forward coefficients {4, 2, 1, 2, 4}.   

--------------------------TASK 03------------------------Check if the system with following filter coefficients is stable.  bk = {1, 3, 5, 3, 2, 2}.   

--------------------------TASK 04------------------------Make a program that checks whether a 3‐point causal averaging filter stable? Plot the pole‐zero  plane. Also determine the frequency response of the filter.

Lab # 13

OBJECTIVES OF THE LAB  

‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐  This lab aims at the understanding of:   • •

Implementation of IIR Filters  Frequency Response of IIR Filters 

‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐‐       

13.1 IIR   IIR stands for infinite impulse response; an example of such a response is step function. In FIR filters, the  output depends on present and previous inputs only. In IIR filters, the output depends not only on  present and previous inputs but also on the previous outputs. We can call FIR filters as special case of IIR  filters. The output is (fed back) to be combined with inputs, hence these systems are also called feed  back systems and the term recursive filters is also used. The z‐transform of IIR filters are rational  functions. They have both poles and zeros at non‐zero locations in the z‐plane     

N

M

l =0

k =0

y[n] = ∑ a l y[n − l ] + ∑ b k x[n − k ]

The filter coefficients consist of two sets of coefficients i.e.   { al } =  coefficients  { bk } = feed forward coefficients   Consider a First order IIR Filter with Filter Coefficients:   

 

 y[n] =0.8 y[n‐1] + 2 x[n] + 2x[n‐1] 

Rewriting:   

 

 y[n] ‐0.8 y[n‐1] = 2 x[n] + 2x[n‐1] 

Therefore,    b = [2, 2]; a = [1, ‐0.8] 

  Example: Frequency Response of the First Order IIR Filter  w = ‐pi:pi/100:pi;  b = [2, 2];  a = [1, ‐0.8];  H = freqz(b,a,w); 

figure;  plot(w, abs(H));  title('MAGNITUDE OF FREQ. RESP. OF FIRST ORDER IIR FILTER y[n] =0.8 y[n‐1] + 2 x[n] + 2x[n‐1]');  xlabel('Normalized Radian Frequency');  ylabel('Magnitude');    figure;  plot(w, angle(H));  title('PHASE OF FREQ. RESP. OF FIRST ORDER IIR FILTER y[n] =0.8 y[n‐1] + 2 x[n] + 2x[n‐1]');  xlabel('Normalized Radian Frequency');  ylabel('Phase (Radians)'); 

 

 

Program: Output of 1st Order IIR Filter with B1=0  %y[n]=0.8y[n‐1] + 5x[n]  b = [5];  a = [1,‐0.8];    % INPUT SIGNAL WITH COEFFICIENTS  x = [2 ‐3 0 2 zeros(1,10)];  y = filter(b,a,x);  s = size(y);  columns = s(2);  n = 0:(columns‐1)  stem(n,y,'linewidth',2.5);  title('OUTPUT SIGNAL FROM A RECURSIVE DIFFERENCE EQUATION y[n]=0.8y[n‐1] + 5x[n]');  xlabel('Discrete‐time variable n');  ylabel('y[n]'); 

-------------------------TASK 01-------------------------Understand how the impulse response of an IIR system is infinite. Take the above mentioned IIR  filter i.e.  y[n]=0.8y[n‐1] + 5x[n] 

Compute by hand (on paper) the values of y for the given input i.e.   

x[n] = 2δ [n] − 3δ [n − 1] + 2δ [n − 3]

And compare it with the results shown by above figure.  Program: Output of 1st Order IIR Filter with B1=0  %y[n] = 1.2y[n‐1] + 5x[n]  b = [5];  a = [1,‐1.2];    % INPUT SIGNAL WITH COEFFICIENTS  x = [2 ‐3 0 2 zeros(1,10)]; 

y = filter(b,a,x);  s = size(y);  columns = s(2);  n = 0:(columns‐1)  stem(n,y,'linewidth',2.5);  title('OUTPUT SIGNAL FROM A RECURSIVE DIFFERENCE EQUATION y[n]=0.8y[n‐1] + 5x[n]');  xlabel('Discrete‐time variable n');  ylabel('y[n]');

Program: Output of 1st Order IIR Filter with B1=0  %y[n] = 1y[n‐1] + 5x[n]  b = [5];  a = [1,‐1];    % INPUT SIGNAL WITH COEFFICIENTS  x = [2 ‐3 0 2 zeros(1,10)]; 

y = filter(b,a,x);  s = size(y);  columns = s(2);  n = 0:(columns‐1)  stem(n,y,'linewidth',2.5);  title('OUTPUT SIGNAL FROM A RECURSIVE DIFFERENCE EQUATION y[n]=0.8y[n‐1] + 5x[n]');  xlabel('Discrete‐time variable n');  ylabel('y[n]');

Program: Output of 1st Order IIR Filter with B1=0  %y[n] = ‐0.5y[n‐1] + ‐4 x[n] + 5x[n‐1]  b = [‐4, 5];  a = [1, 0.5];    % INPUT SIGNAL WITH COEFFICIENTS  x = [2 ‐3 0 2 zeros(1,10)]; 

y = filter(b,a,x);  s = size(y);  columns = s(2);  n = 0:(columns‐1)  stem(n,y,'linewidth',2.5);  title('OUTPUT SIGNAL FROM A RECURSIVE DIFF. EQUATION y[n] =‐0.5y[n‐1]‐4 x[n]+5x[n‐1]');  xlabel('Discrete‐time variable n');  ylabel('y[n]');

-------------------------TASK 02-------------------------Determine the output of an IIR filter when a signal   x = cos(2*pi*0.1*pi*t);  is ed through it. The filter has coefficients  b = [1 2 3 4 5 4 3 2 1];  a = [1 ‐0.1]   

-------------------------TASK 03-------------------------An IIR filter has following coefficients:   bk = {1, 3, ‐2, 4, 1}   al = {‐3, 2, 1}  Input signal   

x[n]= δ[n] ‐  3δ[n‐1] + 3δ[n‐2] + 4δ[n‐3] ‐ 2δ[n‐4]  

 

is applied. Find the output y[n]. 

-------------------------TASK 04-------------------------A recursive discrete‐time system has filter coefficients (feed‐forward) as {1, ‐2, 3, 2} and feed‐ back coefficients as {‐3, 1, 0.5}. An input signal (shown in Fig) is applied to the system.     2 2   1 1 1        ‐2  ‐1    0   1   2   3     4     5     6    7     Determine the discrete‐time output of the filtering. Write code to plot the original signal as well  as the output on two subplots of the same figure.   

--------------------------TASK 05------------------------Plot the magnitude and phase of the frequency response of an IIR filter with coefficients:  b = {1 ‐2 4 ‐2 1} and a = {2, 1, 1, 0.5}