Follow all instructions in the attachments one contains the assginment and the other attachment contains the content and question for the assignment. Detail matlab program is needed.

ECE 115, Sp/13

MATLAB HW #3

Read Chapter 4, titled, “Program Flow Control”, which is posted on the blackboard. Write a

MATLAB program for problem 23 at the end of Chapter 4.

Your program must contain comments stating purpose of the program and explanations of

different parts of the program. Write a well organized and commented program.

At the start of the program, the program must use disp statements to state the program

purpose and how it can be terminated. Then it must get input from the user.

Provide a program listing, and a print out of problem (program) results for several different

inputs.

This is due in class Wednesday, March 6, 2013.

ChapterFOUR

Roland Priemer, Copyright 2012, all rights reserved

___________________________________________________________________

Program Flow Control

A method to solve a problem may include alternative steps. Which next step to take may

depend on results from previous steps. A script that implements a solution method must be

able to follow alternative steps based on answers to questions like: is a variable value different

from an allowed set of values, is the result of a calculation zero, do two variables have the same

value, which of many variables has a particular value, and other questions. The answers to

these kinds of questions can influence what to do next in a script.

After you have completed this chapter, you will know how to:

• test conditions and execute alternative script segments

• repeatedly execute a script segment until a condition occurs

• execute one of many script segments depending on input data or intermediate results.

4.1 Relational Operators

Relational operators are used to compare variable values and produce a result that is true or

false. For example,

>> a = 2; b = 3; % assign values to two variables

>> c = a < b % compare a and b using the relational operator <, meaning less than
c = 1
Since the value of a is less than the value of b, the comparison, a**
**

Table 4.1 List of Relational Operators

Relational Operator Description

< less than

<= less than or equal to > greater than

>= greater than or equal to

== equal to

~= not equal to

Table 4.1 gives the relational operators that can be used in MATLAB. Notice that the

relational operator, ==, which tests equality, is not an assignment operator. The syntax of a

relational expression is

left_expression relational_operator right_expression

The expressions on each side of a relational operator can be any MATLAB expressions that

evaluate to scalars, vectors or matrices having the same dimension. A relational expression

produces a logical result, which can be assigned to a variable.

Example 4.1 ______________________________________________________

>> a = 2; b = 3; c = -1; d = 4; % define four arithmetic variables

>> e = sqrt(a^2+b^2) <= sqrt(c^2+d^2) % arithmetic expressions on each side of the relational operator
e = 1
>> f = (a > b) +2* (c < d) % logical values (and variables) can be used in arithmetic expressions
f = 2
>> % f, an arithmetic variable, can be 0, 1, 2 or 3

>> w = 2*pi; t = 0:0.1:1; % assign a frequency and a time range

>> X = sin(w*t) <= cos(w*t) % compare two vectors element-by-element, X is a logical vector
X = 1 1 0 0 0 0 0 1 1 1 1
>> % matrices that have the same dimensions can be used on each side of a relational operator

>> g = exp(c*t) .* sin(w*t); % define an arithmetic vector, an exponentially decaying sinusoid

>> Z = g > 0 % check times at which exponentially decaying sinusoid is positive

Z = 0 1 1 1 1 1 0 0 0 0 0

% Z is a logical vector, resulting from comparing each element of g to zero

>> g_pos = Z .* g % get positive values of g, a logical vector can be used in arithmetic expressions

g =

Columns 1 through 11

0 0.5319 0.7787 0.7046 0.3940 0.0000 0 0 0 0 0

_____________________________________________________________

Suppose the variables a and b are real scalars. Then, one of the relations

a < b a == b a > b

must be true. However, if a and b are replaced by A and B, which are matrices with the same

dimensions, then the relation is applied to corresponding elements of A and B, and the result is

a logical matrix with the dimension of A, where some elements are logic 1, while others are

logic 0. Since A and B are matrices, then these relations do not produce a scalar logical value.

To test the equality of two matrices, use the built-in MATLAB function isequal, as in c =

isequal(A,B), where c is a scalar logical variable. Table 4.2 gives some built-in MATLAB functions

that operate on matrices and return a logical value.

Table 4.2 Functions That Return a Logical Value

Function Example, A and B are matrices, c is a logical scalar and C is a logical matrix

isequal c = isequal(A,B), are all corresponding elements equal

islogical c = islogical(A), are all elements either 0 or 1

isinteger c = isinteger(A), are all elements integers

isfloat c = isfloat(A), are all elements floating point numbers

isempty c = isempty(A), does A have no elements

true C = true(M,N), C is an MxN logical matrix of all logic 1 values

false C = false(M,N), C is an MxN logical matrix of all logic 0 values

logical C = logical(A), convert numeric matrix A to logical matrix C

any c = any(A), is any element of A nonzero or logical 1, if A is a vector

C = any(A), if A is a matrix, then the function, any, works on each column of A,

and C is a vector. To obtain a scalar, use c = any(any(A))

all c = all(A), are all elements of A nonzero, if A is a vector

C = all(A), if A is a matrix, then the function, all, works on each column of A,

and C is a vector. To obtain a scalar, use c = all(all(A))

Relational expressions evaluate to a logical value. Relational expressions can be

compounded with the logical operators and, or, xor and not.

4.2 Logical Operators

Logical operators work on logical variables and produce a logical result. For example, the

expression, (a < b) and (c == d), can only be true if both relations are true. It is false if either
relation or both relations are false. In MATLAB the word, "and", is replaced by the symbol &,
and a MATLAB assignment statement becomes e = (a**
**

MATLAB uses the symbols &, |, and ~ for the logical operators and, or and not,

respectively, and a built-in function xor for the logical operator xor (exclusive or). Table 4.3

defines the logical operators, where a, b, and c are scalar logical variables.

Table 4.3 Definition of Logical Operators

AND OR XOR NOT

a b c = a & b a b c = a | b a b c = xor(a,b) a ~a

0 0 0 0 0 0 0 0 0 0 1

0 1 0 0 1 1 0 1 1

1 0 0 1 0 1 1 0 1 1 0

1 1 1 1 1 1 1 1 0

Logical operators can work with numbers. If a number is 0, then it is used as if it is

logical 0, and if it is nonzero, then it is used as if it is logical 1.

Example 4.2 ______________________________________________________

>> a = 0; b = 1; c = -3; % assign scalar values

>> a | b % both operands have logical values

ans = 1

>> b & c % c does not have a logical value, but it is used as if it is logic 1

ans = 1

>> ~c

ans = 0

>> xor(a,c) % produces a logic 1 only if one of the two inputs is logic 1

ans = 1

______________________________________________________________

MATLAB also has built-in logical functions for the logical operators. i.e., and(a,b) = a&b,

or(a,b) = a|b, and not(a) = ~a. There is an option when using the & and | operators. Consider

c=expression_a & expression_b. If expression_a evaluates to logic 0 then there is no need to

evaluate expression_b, since c=0. To require that MATLAB evaluates expression_b only if

expression_a is logic 1, use && instead of &. Similarly, consider c=expression_a | expression_b.

If expression_a is logic 1, then there is no need to evaluate expression_b, since c=1. To require

that MATLAB evaluates expression_b only if expression_a is logic 0, use || instead of |. These

options can reduce execution time.

Logical operators can work with matrices when both operands have the same

dimension. The logical operation is done element by element, and the result is a logical matrix.

If one operand is a scalar, then the elements of the other operand, a matrix, are used element

by element to produce a logical matrix. If A is a matrix of numbers, then ~A produces a logical

matrix where a nonzero element of A results in a logic 0 element in ~A. And, ~(~A) produces a

logical matrix, where a nonzero element of A results in a logic 1 element of ~(~A). The built-in

functions given in Table 4.2 can be applied to logical matrices.

Another useful built-in function is the function find to locate the elements of a matrix

that are nonzero. The function returns the indices numbered linearly down each column

starting at the left column.

Example 4.3 ______________________________________________________

>> A = [1 -5 0; -5 0 4; 0 0 2] % assign a matrix

A =

1 -5 0

-5 0 4

0 0 2

>> % the indices for the first column are: 1, 2 and 3, and for the second column

>> % the indices are: 4, 5 and 6, and so on

>> find(A) % find the linear indices of all elements in A that are nonzero

ans = 1 2 4 8 9

>> B = [0 -1.8 -4 5.2 0.6 0 1];

>> I = find(B) % find the indices of the nonzero elements of a vector

I = 2 3 4 5 7

>> find(B,3,’first’) % find the first 3 indices of the nonzero elements in a vector

ans = 2 3 4

>> find(B,4,’last’) % find the last 4 indices of the nonzero elements in a vector

ans = 3 4 5 7

% logical operations and relational operations can be combined in compound expressions

>> find((B > -1) & (B < 1)) % find the indices of the elements in the range: -1 < element < 1
ans = 1 5 6
% Here, (B > -1) produces a logical vector, which is anded element by element with the logic vector

% produced by (B<1) to produce the logical vector that is the input to the find function
______________________________________________________________

Use parentheses to control the sequence in which logical operations are evaluated.

MATLAB evaluates an expression within the innermost parentheses first.

Example 4.4 ______________________________________________________

>> a = 0; b = 1; c = 1; d = 0; % assign logical values

>> a & b | c & ~d % MATLAB evaluates the logical and operations first

ans = 1

% this results from: (a&b) | (c & ~d)

>> a & (b | c) & ~d % MATLAB evaluates the term within parentheses first, and then left to right

ans = 0

% parentheses must be carefully used to achieve the desired overall logical operation

___________________________________________________________________

4.3 If – Elseif – Else – End

With the if-elseif-else structure, a block of statements is executed only if a condition is met.

The syntax of the most general form of this structure is

if expression_1

block #1 of statements

elseif expression_2

block #2 of statements

insert as many additional elseif and block of statements as needed

else

block #N of statements

end

where each expression must evaluate to a logical value. When MATLAB encounters the

if key word, it does the following

1) evaluate expression_1

2) if expression_1 is true, then execute block #1 statements, and continue execution after the

end statement

3) if expression_1 is not true, then skip the block #1 statements to the first elseif key word and

evaluate expression_2

4) if expression_2 is true, then execute block #2 statements, and continue execution after the

end statement

5) if expression_2 is not true, then skip the block #2 statements to the next elseif key word, if

there is one, evaluate its expression, and so on

6) if all expressions are false, then MATLAB will skip to the else key word, if there is one, and

execute the block #N statements, after which execution continues after the end statement.

It is useful to indent each block of statements to easily see that each if key word is matched by

an end key word. The elseif through the else sections are optional. Within each block of

statements there can be nested if-elseif-else-end structures. You can nest these structures as

deeply as needed. The simplest if structure is

if expression

block of statements

end

Here, if the expression is true, then the block of statements is executed, and if it is false, then

the block of statements is skipped and execution continues after the end statement. If the

block of statements can fit on one line, then the if statement could be, for example

>> x = -6.8;

>> % limit x to the range: -5 to +5

if (x<-5) | (x>5), x = 5*sign(x); end % the built-in function sign returns -1, 0 or +1

>> x

x = -5

Example 4.5 _________________________________________________________

Suppose we are given a data point x , where its value can be anywhere in the range

from a to b, written as [ , ]a b . The brackets indicate that the end point values are included in

this range. If we write ( , )a b , then both end points are not included in this range. It is desired

to know which of N sub-intervals of the range contains x . Let us segment this range into N

segments (sub-intervals). Each segment is called a bin. The width of each bin is

( ) /w b a N= − . If we know which bin contains x , then we can say that x is known to be

within its bin range. The centers of the bins are given by

( ) (2 1) , 1, ,

2

w

c k a k k N= + − = (4.1)

Suppose that all we know about x is that it is contained within a particular bin, say bin k .

Then, given k only, the best estimate of the value of x is ( )c k , and the maximum estimation

error is / 2w .

Instead of using bins that all have the same width, a bin center vector c with N

elements could be given, where the bin ranges are: [ , ( (1) (2)) / 2]a c c+ ,

[( (1) (2)) / 2, ( (2) (3)) / 2]c c c c+ + , , [( ( 1) ( )) / 2, ]c N c N b− + , for a total of N bins over [a,b]

that have varying widths. Using bins that all have the same width, Prog. 4.1 finds 1, ,k N= ,

the bin number, given x .

There is another purpose of this script. It illustrates writing a script that is reasonably

robust. This means that input data or intermediate results within the script do not cause a

computer to become non-responsive or generate error messages. A script may be robust and

not give correct (or expected) results, because its programming is based on an incorrect

algorithm. If a script executes to completion and the results are incorrect, then the

programmer (possibly working with a user), must interpret the results, and make corrections.

To find a bin number, the range [ , ]a b is normalized with

x

a

y N

b a

−

=

−

Therefore, if x a= , then 0y = , and if x b= , then y N= , and the range of y is [0, ]N . For a

given x , ceil( )k y= , where ceil is a built-in function that returns the next higher integer.

% Program to find the bin number

clear all; clc

disp(‘Program to find the bin number of x in the range [a,b]’)

err = 0; % used for checking if an error occurred

a = input(‘enter lower range limit: ‘);

b = input(‘enter upper range limit: ‘);

if a < b % a robust program must check inputs
N = input('enter the number of bins, greater than 0: '); N=round(N);
if N > 0 % check input

x = input(‘enter x: ‘);

if x >= a && x <= b % check input
w = (b-a)/N; % get bin width
k = ceil((x-a)/w); % round to next higher integer
else
err=1; % x is out of range
end
else
err=2; % N is out of range
end
else
err=3; % lower limit not less than upper limit
end
if ~err % was an error detected
fprintf('out of N = %i bins, the bin number is: %i \n',N,k)
elseif err == 1, disp('x is out of range')
elseif err == 2, disp('N is out of range')
elseif err == 3, disp('lower limit not less than upper limit')
end

Program 4.1. Program to find the bin number of x in the range [a,b].

Notice that all statements within matched if, else, elseif, and end key words are indented for

clarity.

___________________________________________________________________

4.4 For Loop

With the for loop structure, you can repeat executing a block of statements a predetermined

number of times. Its syntax is

for loop_index = first_index : index_increment : last_index

block of statements

end

Every for key word must be matched by an end key word. The block of statements can include

additional for loops. For loops can be nested. The block of statements can include if-elseif-

else-end structures and any other valid MATLAB statements. It is useful to indent the block of

statements to easily see the matched for and end key words. The loop_index must be a

MATLAB variable. The first_index value, index_increment value and last_index value can be

expressions that evaluate to a real scalar.

For convenience, denote the loop_index variable by k, first_index by a, index_increment

by incr, and last_index by b. When MATLAB encounters the for key word it does the following

1) check incr, and if incr is positive, check that a <= b, and if a <= b, then set k = a 2) execute the block of statements 3) go back to the for statement and set k = k + incr, which gives the new loop_index value 4) if k <= b, then execute the block of statements again and go back to step (3) 5) if k > b, then skip the block of statements, and continue execution after the end statement.

If a>b, then the block of statements is not executed, and execution continues after the end

statement.

If in step (1), incr is negative, then check that a >= b, and if a >= b, then set k = a, and do

the following

2) execute the block of statements

3) go back to the for statement and set k = k + incr, which gives the new loop_index value

4) if k >= b, then execute the block of statements again and go back to step (3)

5) if k**
**

where v is a vector and k is first set to v(1) and the block of statements is executed . Then,

k=v(2) and the block of statements is executed, and this continues until all elements of v have

been used. With this structure you can assign to k any values, including complex values, for

example, v = [1-j*2 2* pi -1 exp(j*pi/4) u], where u = [1/3 j*3].

If you do not include an index_increment, as in

for k = a:b

block of statements

end

then a default index_increment value of 1 is used, which requires that a<=b. Example 4.6 _______________________________________________________ Let us take a look at the details of multiplying two matrices with for loops. Let A be an (NxK) matrix, and let B be a (KxM) matrix. Consider C = A*B with dimension (NxM). Prog. 4.2 finds C with for loops and the built-in MATLAB matrix multiply operation. % Program to multiply two matrices (C=A*B) using for loops clear all; clc N=3; K=4; M=2; % specify matrix dimensions A=[1 0 2 -1;2 1 -1 0; 1 -1 1 2]; % specify two matrices B=[3 0; -1 5; 1 3; 2 1]; C = zeros(N,M); % preallocating space for C, which decreases execution time t = tic; % tic, a built-in function, starts a stopwatch timer for n = 1:N % index for rows of C for m = 1:M % index for columns of C sum = 0; % initializing the sum of products for k = 1:K % get dot product of a row of A and a column of B sum = sum + A(n,k)*B(k,m); % accumulate products end C(n,m) = sum; % save dot product of a row of A and a column of B in C end end loops_compute_time = toc(t) % toc gets elapsed computing time since tic C % display result t = tic % start stopwatch timer D=A*B; % check MATLAB_compute_time = toc(t) % get elapsed computing time since tic

Program 4.2. Use of for loops to multiply two matrices.

Notice that all statements within matched for and end key words are indented for

clarity. Since the innermost loop has only a few statements, it could have been written as

for k = 1:K, sum = sum + A(n,k)*B(k,m); end

This for loop could also have been replaced to find C(n,m) with

C(n,m) = A(n,:)*B(:,m); % dot product of row n in A and column m in B

The script output is shown below. It was rearranged to take less space.

loops_compute_time = 0.0024 secs

C =

3 5

4 2

9 0

MATLAB_compute_time = 1.4848e-005 secs

This example demonstrates a programming utility of MATLAB, i.e., it works with the

matrix data type. Moreover, D=A*B is executed in much less computing time than it took to

obtain C, even though space for C was preallocated.

___________________________________________________________________

The built-in functions tic and toc were used in Example 4.6 to measure the execution

times of parts of Prog. 4.2. The time required to execute a part of a script depends on the

speed of your computer. There may be occasions when, for example, between different

program outputs, you want to control the time in between these program outputs. To do this,

use the built-in function pause(seconds), for example, pause(4.5), which causes a program to

stop execution for 4.5 secs. Use pause without an argument to stop a program and wait for the

user to strike any key before continuing.

If in a script the size of a matrix (or vector) is increased with each pass through a loop, it

is best to preallocate space for the largest size of the matrix before getting into the loop. This

will save much execution time.

Example 4.7 ______________________________________________________

A signal is ( ) sin( )x t A tω= , and we want to see how the values of ( )x t are distributed

over one cycle. The signal value varies over the range A− to A , written as [ , ]A A− . Let us

segment this range into N bins. The width of each bin is 2 /A N , and using (4.1) the bin

centers are given by

( ) (2 1) , 1, 2, ,

A

c k A k k N

N

= − + − = (4.2)

Let 2ω π= rad/sec, and therefore one cycle takes 0 1T = sec. The signal ( )x t will be sampled at the

rate 44,100sf = samples/sec. Prog. 4.3 counts the number of times that the value of ( )x t is in

each of the N bins. The resulting distribution is called a histogram.

% Program to find the histogram of a sine wave

clear all; clc

% specify the data

A = 1; freq = 1; T0 = 1/freq; % specify amplitude and freuency

fs = 44100; T = 1/fs; % sampling rate and time increment

t = 0.0:T:T0-T; % time points over one cycle

x = A*sin(2*pi*freq*t); N_data = length(x); % get samples of x(t)

N = 40; w = 2*A/N; k = 1:N; % number of bins, bin width, and bin index

c = -A+(2*k-1)*w/2; % bin centers

bin = zeros(1,N); % preallocate space and initialize bin counts to zero

for n=1:N_data

k = ceil((x(n)+1)/w);

if k==0

k=1; % x=-1 belongs to bin 1

end

bin(k)=bin(k)+1; % increment the bin count

end

bar(c,bin,1) % get a barchart plot

grid on

% the label were entered using the plot editor

Program 4.3. Program to find the histogram of a sine wave.

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8

1

0

500

1000

1500

2000

2500

3000

3500

4000

4500

sine wave value

va

lu

e

di

s

t

rib

ut

io

n

sine wave histogram

Figure 4.1. Histogram of a sinusoidal signal.

The histogram is shown in Fig. 4.1. We see that a sine wave spends almost equal amounts of

time in segments from -0.8 to +0.8, while nearly 20% of the time in a cycle is spent within 5% of

the maximum and minimum values.

This example introduced the built-in function bar, which plots a bar chart. Like the built-

in function plot, the first input gives horizontal plotting information, the bar centers, the second

input gives the bar heights, and the third input controls spacing between bars in the plot, where

a value less than 1 causes a space between bars. More detail will be given in Chapter 9.

______________________________________________________________

In the previous examples, the for loops always executed to completion. Sometimes it is

useful to terminate a for loop before it has iterated through its entire index range.

Example 4.8 ______________________________________________________

Evaluate the sin( )x function using its Taylor series, which is given by

2 1

0

sin( )

( 1)

(2 1)!

k k

k

x

x

k

∞ +

=

=

−

+∑

In a script, the upper summation limit cannot be infinity. We must either set the upper

limit to a large integer or use some way of stopping the computation after some predetermined

precision has been achieved. When x is incremented by 2π the function sin( )x repeats itself.

Let us evaluate the function for 0 2x π≤ ≤ . Program 4.4 evaluates the Taylor series and plots

the result. The resulting approximation of the sin function is shown in Fig. 4.2.

% Program to evaluate the sin(x) function using its Taylor series

clear all; clc

N = 101; % N values for x will be used

theta = linspace(0,2*pi,N); % assign values to theta

y = sin(theta); % use the built-in function sin to compare with the series

K = 100; % K terms in the Taylor series will be used

precision = 1e-6; % this is a precision of 0.0001 percent

y_approx = zeros(1,N); k_value = zeros(1,N); % preallocate space

for n=1:N % loop for each value in theta

x = theta(n);

% no need to evaluate the Taylor series for x=0

if x == 0.0, y_approx(n) = 0; k_value(n) = 0; continue; end

sum_T = 0; % initialize Taylor series sum

factorial_term = 1; % initialize (2k+1)!

minus_one_term = 1; % initialize (-1)^k

x_power_term = x; % initialize x^(2k+1)

x_square = x^2; % used to get next x^(2k+1), computed only once

next_term = minus_one_term*x_power_term/factorial_term; % k=0 term

for k=1:K % loop to accumulate series terms

sum_T = sum_T + next_term; % accumulate terms

% prepare the next term in the Taylor series

factorial_term = factorial_term*(2*k)*(2*k+1); % update (2k+1)!

% uses only two multiplies to find each (2k+1)!

minus_one_term = – minus_one_term; % update (-1)^k

x_power_term = x_power_term*x_square; % update x^(2k+1)

% uses only one multiply to find each x^(2k+1)

next_term = minus_one_term*x_power_term/factorial_term;

% compare next term in Taylor series to sum and check precision

if abs(next_term) >= abs(precision*sum_T);

continue % continue this for loop

else

sum_T = sum_T + next_term; % include the next term

break % jump out of this for loop

end

end

% the size of y_approx is increased with each pass through this loop

% space for its largest size was preallocated before this loop

y_approx(n) = sum_T;

k_value(n) = k; % save index for which this for loop achieved precision

end

% finished evaluating the Taylor series for each element of theta

max_k_value = max(k_value) % get maximum number of for loop iterations

max_error = max(abs(y-y_approx)) % get maximum error

plot(theta,y_approx)

grid on

xlabel(‘angle – radians’)

title(‘sin function with Taylor series’)

Program 4.4. Evaluate the sin function using its Taylor series.

0 1 2 3 4 5 6 7

-1

–

0.

5

0

0.5

1

angle – radians

sin function with Taylor series

Figure 4.2. sin function produced with its Taylor series.

The outer for loop in Prog. 4.4 iterates over the values in the vector theta. The inner for loop in

Prog. 4.4 is set up to iterate K=100 times. It computes each term in the Taylor series. Since

K=100 terms may not be necessary to achieve a desired precision, each new term is checked to

see if it has become small enough to make a negligible contribution to the summation. If it has

not become small enough, the built-in function continue is used to continue the inner for loop

iteration. If it has become small enough, then the built-in function break is used to exit the

inner for loop in which break occurs. Some script parameters of interest are shown below.

While the inner for loop could have iterated up to K=100 times, the maximum number of

iterations to achieve the prescribed precision of 0.0001 % over a cycle of the sin function

turned out to be k=22, and the resulting maximum error when compared to sin( )x using the

MATLAB built-in function is 3.7e-6 percent. The use of tic and toc will show that the algorithm

that MATLAB uses to evaluate the sin function is much faster than a direct application of its

Taylor series.

max_k_value = 22

max_error = 3.7065e-008

______________________________________________________________

Some built-in functions concerned with controlling loop activity are given in Table 4.4.

Table 4.4 MATLAB Functions to Exit For Loops, While Loops and User-Defined Functions

Function Brief Description

continue When continue is encountered in a loop, remaining loop statements are

skipped and execution continues with the next iteration of the loop in which

it appears. In nested loops, continue causes the next iteration of the loop

enclosing it.

break When break is encountered in a loop, it terminates execution of the loop. In

nested loops, break cause an exit from the loop enclosing it only.

return When return is encountered in a user-defined function, it causes an exit from

the function. Statements in the function after the return statement are not

executed. Also, terminates the keyboard mode, which will be discussed

Chapter 10.

While loops will be discussed in the next section.

4.4.1 Probability

To consider some interesting problems, a brief background in probability will now be given.

Suppose we do an experiment, and the outcome is not predictable. Assume that an

outcome can be one particular event among a finite number of possible events. The probability

of an outcome gives an assessment of how likely it is that a particular outcome will occur

compared to all other possible outcomes.

For example, consider tossing a coin. Each coin toss is an experiment, and there are two

possible outcomes or events, i.e., heads or tails. Suppose you toss a coin N times. Of the

N tosses, let Hn be the count of the outcome heads, and let Tn be the count of the outcome

tails. The probability Hp of heads is given by

lim HH N

n

p

N→ ∞

= (4.3)

and for tails, the probability Tp is given by

lim TT N

n

p

N→ ∞

= (4.4)

Since H TN n n= + , we have a fundamental property of the probabilities of all possible

outcomes, which in this case is:

1lim limH T H TN Np p

n n

N N→∞ →∞

+ = +

=

For a fair coin we expect that 0.5Hp = and 0.5Tp = .

If there are K possible outcomes of an experiment, and , 1, 2, ,kp k K= is the

probability of each outcome, then

1

1

K

k

k

p

=

=∑ (4.5)

Practically, we cannot do an experiment an infinite number of times to find a probability

by counting. However, if we make N large enough, then we might get useful values for the

probabilities. The probability of an impossible outcome is zero, while the probability of a

certain outcome is one. For the probability p of any outcome we have 0 1p≤ ≤ . Notice that if

an outcome occurs only once in an infinite number of trials of an experiment, then its

probability is zero. This means that if the probability of an outcome is zero, we cannot be

certain that it will never occur. Similarly, if the probability of an outcome is one, then we

cannot be certain that it did not occur at least once.

We can use the MATLAB built in function rand to experimentally find probabilities

associated with the outcomes of some experiments. In MATLAB, each rand function evaluation

returns a number in the range: 0 1→ (written as (0,1)), which we cannot predict, and where

all numbers in this range are equally likely, which means that the function rand is a uniformly

distributed random number generator. The MATLAB statement

x = rand; % 0 < x < 1 assigns to x a number in the range (0,1). Each use of rand returns a number that does not depend on numbers obtained previously. A pseudorandom number generator is an algorithm for generating a sequence of numbers that approximates the properties of random numbers. The sequence is not truly random. The function rand is a pseudorandom number generator. A careful mathematical analysis is required to have any confidence that a pseudorandom number generator produces numbers that are sufficiently random to meet the needs of a particular application. Nevertheless, it is used in simulations to study the behavior of systems in which random events can occur.

http://en.wikipedia.org/wiki/Algorithm�

http://en.wikipedia.org/wiki/Random�

The algorithm for the function rand is initialized the same way each time a MATLAB

session is started. Therefore, when you start a MATLAB session, rand will produce the same

random sequence of numbers that it produced the previous time that a MATLAB session was

started. During a MATLAB session you can reset the random number generator to the default

MATLAB start-up state with the statement: rng(‘default’). Use the MATLAB help facility to find

out more about the rand function.

Example 4.9 _____________________________________________________

Let us simulate the coin toss experiment to find Hp and Tp . This will be done by

counting the number of times each outcome occurs over a large number of coin tosses. Prog.

4.5 simulates the coin toss experiment.

clear all; clc;

rng(‘default’); % initialize random number generator to its default state

total_tosses = 1e6; % specify the total number of coin tosses

heads = 0; tails = 0; % initialize outcome counters

for n = 1:total_tosses

% do experiment

x = rand; % get a random number distributed from 0 to 1

if x < 0.5 % let this result be analogous to heads occurring
heads = heads + 1;
else % if not heads, then tails
tails = tails + 1;
end
end
disp('probability of heads: ');
pH = heads/total_tosses
disp('probability of tails: ');
pT = tails/total_tosses

Program 4.5. Using a random number generator to simulate tossing a coin.

The program output is shown below.

probability of heads: pH = 0.4995

probability of tails: pT = 0.5005

If we increase the number of tosses, then it is likely that pH and pT will each be closer to 0.5.

However, this does not mean that as the total number of tosses N is increased, then the

number Hn of heads or the number Tn of tails will ever become exactly / 2N . Since pH is

essentially the same as pT, the outcomes heads and tails are equally likely.

_______________________________________________________________

The mean (average) value of the numbers produced by rand can be found

experimentally with

rng(‘default’) % initialize random number generator

N = 1e7; % specify the number of random numbers to get

% obtain a vector of N random numbers uniformly distributed from 0 to 1

x = rand(1,N);

x_mean = sum(x)/N % find the mean value

% this can be obtained with the MATLAB function mean to get x_mean=mean(x)

x_mean = 0.5000

With the function rand we can obtain uniformly distributed random numbers over any

range (a,b). The span of this range is (b-a). Since the range of x = rand is (0,1), the range and

span of y = (b-a)x are (0,(b-a)) and (b-a), respectively. Therefore, the range of y = (b-a)x + a is

(a,b).

Example 4.10 _____________________________________________________

To obtain random numbers x with a uniform distribution over the range: -5 < x < +5, use
>> N = 1e7; % specify the number of random numbers

>> x = rand(1,N); % get a vector of N uniformly distributed random numbers between 0 and 1

>> x = 10*x – 5; % multiply by (5-(-5))=10 and translate by -5

>> x_mean = mean(x) % find the mean value

x_mean = 2.9411e-006

>> % ideally, the mean value of x should be x_mean = 0.

>>

Let us see the first 100 points in x with

>> y = x(1:100); % get first 100 points from x

>> y = y – mean(y); % make the mean of y zero

>> plot(y,’- .’); grid on; xlabel(‘index’); ylabel(‘random signal’)

>> % within single quotes, the dash causes the plot function to connect the points with straight

>> % lines, and the period causes the plot function to put a dot at each point

More about plotting will be given in Chapter 9. Fig. 4.3 shows a plot of some points from x. It

looks like a noise signal with a zero mean value.

0 10 20 30 40 50 60 70 80 90 100

-5

0

5

index

ra

nd

om

s

ig

na

l

Figure 4.3. A random signal.

______________________________________________________________

An important communications problem occurs due to noise. This is because a received

signal ( )x t is a noise corrupted version of a transmitted signal ( )s t , such as

( ) ( ) ( )x t s t n t= +

where ( )n t is additive noise. For example, the signal ( )x t could be the playback from an old

audio recording, and ( )s t is the audio signal that was originally recorded. Or, ( )x t could be a

received cell phone signal, and ( )s t is the signal that was transmitted from a distant location.

We are interested to design a process that performs as depicted in Fig. 4.4, where the process

output ˆ( )s t is an estimate of ( )s t .

( )s t

( )x t

Process that operates on

its input x(t) to produce an

estimate of s(t)

Figure 4.4. Pictorial of a desired operation.

To design an estimation process, it will be useful to do simulation studies with a computer to

test the performance of various process designs and noise signals. This will require simulating

the noise signal with a model of noise, such as the noise signal shown in Fig. 4.3. Such

simulations are often done with MATLAB.

MATLAB also includes the two built-in functions randi and randn. The function randi

produces uniformly distributed random integers. For example, A = randi(imax,M,N) is an (MxN)

matrix of random integers, where each element in A is in the range 1 to imax, while x =

randi(imax) is a scalar integer in the range 1 to imax.

Before describing the random numbers produced by randn, let us study the random

numbers it produces by finding the mean, variance and a histogram. The variance, denoted by

2σ , of a set of random numbers is a measure of their volatility. The variance is computed with

2 2

1

1 ( ( ) )

N

n

x n x

N

σ

=

= −∑ (4.6)

where x denotes the mean of x , and σ is called the standard deviation. More specifically,

the variance (and standard deviation) of a set of random numbers is a measure of the degree to

which the random numbers in the set deviate away from their mean. MATLAB has a built-in

function, var, that computes the variance. A histogram of a set of random numbers shows how

the numbers are distributed.

The built-in function hist returns a vector, the elements of which are counts of the

number of times the elements of the vector of random numbers have values within the range

of each of the bins. We could have used the hist function in Example 4.7.

Example 4.11 _____________________________________________________

Prog. 4.6 finds the mean, variance and a histogram of the numbers obtained with

randn.

% Histogram of numbers produced by the psuedo random number generator randn

clear all; clc

N = 1e7; % specify the number of random numbers to get

x = randn(1,N); % obtain a vector of N random numbers

x_mean = mean(x) % find the mean value

x_var = var(x) % get the variance of the random numbers in x

bin_width = 0.25; % specify a bin width less than 1

N_pos_x_bins = 25; % specify number of bins for positive x(n) values

pos_x_bin_center = N_pos_x_bins*bin_width; % right most bin center

% locate equally spaced bin centers

bin_centers = -pos_x_bin_center: bin_width : pos_x_bin_center;

fprintf(‘using %i bins \n’,length(bin_centers))

x_hist = hist(x,bin_centers); % return a vector with the count in each bin

% prob_x = x_hist/N; % probability of x within each bin

bar(bin_centers,x_hist,1,’w’); % full bar width; and white color

axis([-4 4 0 1.05*max(x_hist)]); grid on

xlabel(‘random number value’); ylabel(‘bin count’)

title(‘Distribution’);

bins_to_1 = round(1/bin_width); % number of bins accounting for 0

% this is the probability that a value of x will be in the range –sigma to

% +sigma, where sigma is the standard deviation. With randn, sigma = 1

% additional plot editing was done in the Figure Window

Program 4.6. Script to find the mean, variance and a histogram of a set of numbers.

The program output, which is given below, shows that the built-in function randn produces

random numbers having a zero mean and a variance of one. The standard deviation is one, and

within ± one standard deviation 67.92 % of the returned values occur. The histogram has a bell

shaped appearance. A normally (also called Gaussian) distributed random number sequence

has these properties.

x_mean = -5.9235e-004

x_var = 0.9993

using 51 bins

prob_minus_to_plus_sigma = 0.6792

The function randn is designed to be a normally distributed pseudorandom number

generator, and produces numbers with zero mean and unity variance. With y a x b= + you

can produce a normally distributed random number sequence y having a prescribed mean and

variance, where 2 2y aσ = and y b= .

-4 -3 -2 -1 0 1 2 3

4

0

2

4

6

8

10

x 10

5

random number value

bi

n

co

un

t

Distribution

Figure 4.5. Distribution of random numbers obtained with the function randn.

If each bin count is divided by the total of number of samples N used to get the

histogram of Fig. 4.5, then, like (4.3) and (4.4) , we get the probability that a number will occur

in a particular bin. Designate the bin counts for positive (negative) numbers by kb ( kb− ),

respectively. Then, for example, 0b is the bin count for 0.125 0.125x− ≤ < , and the probability
that an experiment outcome (invocation of randn) will result in a value of x within the range

0.125 0.125x− ≤ < is 0 0 / 10 5 / 1 7 0.1p b N e e= = . The function randn produces numbers in

the range x−∞ < <+ ∞ , where a large value of abs( )x is very unlikely, and like (4.5) we have that

1k

k

p

+∞

=−∞

=∑ (4.7)

______________________________________________________________

4.4.2 Median Filtering

In signal processing, it is often desirable to perform some kind of noise reduction on a signal

(one-dimensional or 1-D signal) or an image (2-D signal). Median filtering is a digital signal

processing algorithm that is easy to apply and often used to reduce noise. For example, an old

audio recording, where there is crackling noise, can be improved by median filtering. Median

filtering of a speech signal is sometimes a pre-processing step before speech recognition

algorithms are applied. Median filtering is widely used in digital image processing.

The median filter concept is actually simple. Let us see how it works with a

demonstration. Suppose a signal s(t) is sampled, like in (3.33), to produce the number

sequence (discrete time signal) x=[5 6 47 7 6 -29 4], where x(1)=s(0), x(2)=s(T), x(3)=s(2T),

and so on. The third and sixth elements of x appear to be unusual, and they are called outliers.

Let y denote the number sequence produced by median filtering x. For element n of y, use x(n-

1), x(n) and x(n+1). This is called a K=3 element window, w = [x(n-1) x(n) x(n+1)]. Sort w into

ascending values to obtain a K=3 element vector u. Now set y(n) to the middle element of u.

Then, increment n by 1, get the next w, and repeat to get the next y. For example,

n=1, w=[5 5 6] → sort → u=[5 5 6], y(1)=u(2)=5. There is no x(0); x(1) was repeated.

n=2, w=[5 6 47] → sort → u=[5 6 47], y(2)=6

n=3, w=[6 47 7] → sort → u=[6 7 47], y(3)=7

n=4, w=[47 7 6] → sort → u=[6 7 47], y(4)=7

n=5, w=[7 6 -29] → sort → u=[-29 6 7], y(5)=6

n=6, w=[6 -29 4] → sort → u=[-29 4 6], y(6)=4

n=7, w=[-29 4 4] → sort → u=[-29 4 4], y(7)=4. There is no x(8); x(7) was repeated.

The outliers in x have been removed by median filtering. The signal y is less erratic (noisy) than

the signal x. Median filtering works best to replace outliers with a neighboring element when

other consecutive elements in x have similar values.

The window w is said to slide across the signal x. For this demonstration K=3 was used.

We could have used K=4, in which case w is set to w=[x(n-2) x(n-1) x(n) x(n+1)], and u will not

have a middle element. In this case, set y(n)=(u(n-1)+u(n))/2. Or, we can apply the K=3 median

filter again to y, and the result will be even less noisy and most likely too different from x

without the noise. The built-in MATLAB function median finds the median of a vector w.

http://en.wikipedia.org/wiki/Signal_processing�

http://en.wikipedia.org/wiki/Noise_reduction�

http://en.wikipedia.org/wiki/Digital_filter�

http://en.wikipedia.org/wiki/Digital_filter�

http://en.wikipedia.org/wiki/Signal_noise�

http://en.wikipedia.org/wiki/Image_processing�

Example 4.12 _________________________________________________________

For an 88 key piano, let 1, 2, , 88n = be the key number from the lowest to the

highest frequency tone. The frequency (pitch) f of the tone produced by key n is given by

1/12 ( 49)440 (2 ) nf −= Hz (4.8)

For 40n = the note is middle C, and for 49n = , the note played is A440.

Prog. 4.7 consists of two parts. In the first part, a user can enter an integer n , and then

the program uses (4.8) to get the frequency of the note. A sinusoid at this frequency is then

sampled at the rate of 44,100 samples/sec for 5 secs. To make the sound produced more

realistic, the sampled sinusoid is multiplied by the pulse shown in Fig. 3.3, which becomes the

envelope of the sound. Therefore, the sound increases quickly and decays slowly. To hear the

sound, a two column (two channel stereo) matrix must be set up, and then, the built-in

MATLAB functions sound or soundsc can be used to send the matrix to the audio play device of

the computer. You can also use the built-in MATLAB function wavwrite, which creates a

Windows WAV file that can be played by any multimedia application. The syntax for wavwrite

is wavwrite(x,fs,’name.wav’), where x is a two column matrix, fs is the sampling rate (also called

sampling frequency), and name can be replaced by a name you can choose.

% Program to sound an 88-key piano

% Keys are identified by an integer n, n = 1 to 88 from left to right

% Only the tone at the fundamental frequency (pitch) will be sounded

clc; clear all

fs = 44100; % this is the sampling rate used to produce audio CDs.

disp(‘Play a note on a virtual piano.’)

n = input(‘enter an integer from 1 to 88 of a piano key: ‘); % using n=40

n = round(n); % make sure entry is an integer

if n < 1, n = 1; end % keep n in range
if n > 88, n = 88; end

f = 440*(2^(1/12))^(n-49); % the frequency in Hz of the key

w0 = 2*pi*f; % convert Hz to rad/sec

T0 = 5; % the tone will be sounded for 5 secs

delta_t = 1/fs; % the time increment between sound sample

N = round(T0/delta_t); % integer number of samples

time = delta_t*[0:N-1]; % N time points where sound is sampled

p_key = sin(w0*time); % get samples of sound, limit amplitude to +/-1

% get an envelope

load(‘pulse.mat’,’x’) % retrieve the pulse that was built in Example 3.15

p_key = x .* p_key; % multiply the sound by an envelope, this is optional

stereo_sound = [p_key’ p_key’]; % form left and right channels

sound(stereo_sound,fs) % send sound to audio play device

Program 4.7 (first part). Sound a piano key.

If the note produced in the first part of Prog. 4.7 comes from an old recording, then it

will likely contain some spike like noise. The second part of Prog. 4.7 simulates the addition of

noise and spike noise, which occurs randomly in the sound, and then it applies median filtering

to reduce the spike noise contribution to the sound.

% example of median filtering

noise = 0.1*randn(1,N); % using the normal random number generator

p_key = p_key + noise; % add noise to signal

% generating randomly positioned spikes

increment = 100; % add noise spike somewhere within every 100 sample of p_key

spikes = zeros(1,N); % zero out and preallocate space

for k = 1:increment:N;

location = floor(k + rand*increment); % use rand to locate spike

if rand > 0.5; % use rand to determine sign of the noise spike

sign = +1.0;

else

sign = -1.0;

end

spikes(location) = sign*2.0*rand; % use rand to determine spike value

end

p_key = p_key + spikes; % add spike noise to signal

plot_points = 1500; range = 5e4:5e4+plot_points-1; % plotting few points

t_interval = time(range); y = p_key(range);

plot(t_interval,y);

axis([t_interval(1) t_interval(plot_points) -3 3])

title(‘Segment of note middle C plus random noise plus random spikes’);

xlabel(‘seconds’);

ylabel(‘amplitude’);

grid on

max_value = max(abs(p_key));

% limit sound to +/- 1, required by audio device

normalized_p_key = p_key/max_value;

stereo_sound =[normalized_p_key’, normalized_p_key’];

sound(stereo_sound,fs)

% median filtering sound

filtered_p_key = zeros(1,N); % preallocate space

K = 3; % use 3 element median window

for k = 1:N

if k == 1

w = [p_key(1) p_key(1) p_key(2)];

filtered_p_key(1) = median(w);

elseif k == N

w = [p_key(N-1) p_key(N) p_key(N)];

filtered_p_key(N) = median(w);

else

w = [p_key(k-1) p_key(k) p_key(k+1)];

filtered_p_key(k) = median(w);

end

end

y = filtered_p_key(range); % plot a segment of filtered_p_key sound

figure % causes function plot to open a new figure window

plot(t_interval,y)

axis([t_interval(1) t_interval(plot_points) -3 3])

title(‘Segment of median filtered note middle C’);

xlabel(‘seconds’);

ylabel(‘amplitude’);

grid on

max_value = max(abs(filtered_p_key));

% limit sound to +/- 1, required by audio device

normalized_filtered_p_key = filtered_p_key/max_value;

stereo_sound =[normalized_filtered_p_key’, normalized_filtered_p_key’];

sound(stereo_sound,fs)

Program 4.7(second part). Median filter piano key sound.

Fig. 4.6 shows note middle C with noise and spike noise, and the sound produced by this

signal is awful. Fig. 4.7 shows the result of median filtering to remove spike noise.

1.135 1.14 1.145 1.15 1.155 1.16 1.165

-3

-2

-1

0

1

2

3

Segment of note middle C plus random noise plus random spikes

seconds

am

pl

itu

de

Figure 4.6. Note middle C with noise and spike noise.

1.135 1.14 1.145 1.15 1.155 1.16 1.165

-3

-2

-1

0

1

2

3

Segment of median filtered note middle C

seconds

am

pl

itu

de

Figure 4.7. Median filtered note middle C.

The signal shown in Fig. 4.7 still has some noise in it. However, additional digital signal

processing can be performed to reduce the noise further.

__________________________________________________________________

4.5 While Loop

With the while loop structure you can repeat executing a block of statements until an

expression is no longer true. Its syntax is

while expression

block of statements

end

Every while key word must be matched by an end key word. The block of statements can

include additional while loops and for loops. While loops can be nested. The block of

statements can include if-elseif-else-end structures and any other valid MATLAB statements. It

is useful to indent the block of statements to easily see the matched while and end key words.

The expression must evaluate to a logical or real scalar. If it evaluates to a real number, then it

is treated as a logical value.

When MATLAB encounters the while key word it does the following

1) evaluate the expression

2) if the expression is true, then execute the block of statements

3) and go back to the while statement and continue with step (1)

4) if the expression is not true, then skip the block of statements, and continue execution after

the end statement.

A while loop can be made to work like a for loop. For example,

n = 1; N = 10;

while n <= N
n = n +1;
end
Here, the statements within the while loop will be executed N times. However, this is not the
intention of including a while loop structure in addition to a for loop structure in the
programming language. Usually, a while loop structure is used when we do not know ahead of
time how many loop iterations are necessary to accomplish some task. For example,
err = 0;
while ~err
block of statements

end

Within the block of statements, some condition must be tested that will either leave err

unchanged and the loop continues to iterate or set err to err=1, which causes the while loop to

terminate.

Example 4.13 _________________________________________________________

The inner for loop in Prog. 4.4 of Example 4.8 could be replaced with

k = 1; negligible = 0;

while ~negligible

sum_T = sum_T + next_term; % accumulate terms

% prepare the next term in the Taylor series

factorial_term = factorial_term*(2*k)*(2*k+1); % update (2k+1)!

% uses only two multiplies to find each (2k+1)!

minus_one_term = – minus_one_term; % update (-1)^k

x_power_term = x_power_term*x_square; % update x^(2k+1)

% uses only one multiply to find each x^(2k+1)

next_term = minus_one_term*x_power_term/factorial_term;

% compare next term in Taylor series to sum and check precision

if abs(next_term) < abs(precision*sum_T);
sum_T = sum_T + next_term; % include the next term
negligible = 1; % no need for additional loop iterations
end
k = k + 1; % count iterations required to achieve precision
end
The if-end block of statements checks a condition that can change the result of evaluating the
while expression.
______________________________________________________________
The break and continue commands can also be used in while loops. They work the
same way as in for loops. If a while loop expression is always true, then, for example, we have
while 1
block of statements
end
The block of statements must include a test of some condition that can cause a break,
otherwise the loop will execute endlessly. You can terminate a MATLAB program by activating
the Command Window and depressing Ctrl and C on the keyboard.
4.6 Method of Steepest Descent

Consider the problem of finding the value of x that minimizes a function ( )f x , assuming that

( )f x has a minimum. Denote the value of x that minimizes ( )f x by optx . Therefore,

( ) ( )optf x f x≤ , for optx x≠ . A method to find optx is based on the drawing shown in Fig. 4.8.

We can see three minima in Fig. 4.8. At optx x= , ( )f x has its global minimum (also called the

optimizer of ( )f x ), and the other two minima at ax x= and bx x= are called local minima.

f(x)

xx opt x nx bx a

Figure 4.8. A function with a minimum at x = xopt.

Suppose nx is a guess of optx . In Fig. 4.8, nx is to the right of optx , where the slope of

( )f x at nx x= is positive, as indicated by the tangential line. Therefore, if the slope of ( )f x

at the guess nx is positive, then we know that we must decrease this guess. Similarly, if nx is to

the left (but not too far) of optx , then the slope of ( )f x at nx will be negative, and we know

that we must increase this guess. Therefore, given a guess nx of optx , adjust the value of nx to

obtain the next guess 1nx + with

‘1 ( )n n nx x f xµ+ = − (4.9)

where µ , called the step size, is a small positive number that we must pick and ‘( )nf x is the

derivative of ( )f x , assuming it exists, evaluated at nx x= .

To find a value of x that maximizes a function ( )f x , change (4.9) to

‘1 ( )n n nx x f xµ+ = + (4.10)

The recursion, (4.9), of the method of steepest descent is operated for 0, 1, 2,n = ,

where 0x = initial guess of optx . If 0x is within the vicinity of optx and µ is small enough, then

lim n optn x x→∞ = , since as n optx x→ , ‘( ) 0nf x → . Selecting a suitable value of µ is a matter of

trial and error, and depends on the behavior of ( )f x as well as computation precision. If the

step size is too big, then the next value 1nx + may be further away from optx than nx , or even in

the vicinity of a local minimum, such as ax or bx in Fig. 4.8, in which case (4.9) may never

converge to optx . If the step size is too small, then (4.9) will require many iterations (much

computing time) to converge. The recursion (4.9) (or (4.10)) is operated until ε+ − <1n nx x , for some desired precision ε . We do not have to keep the sequence nx , 0, 1, 2,n = . Equation (4.9) could be an update of oldx to obtain newx , as in 1) set µ to a small positive number and set oldx to an initial guess of optx

2) compute ‘ ( )new old oldx x f xµ= −

3) if ε− ≥new oldx x , then old newx x= and repeat step (2), else opt newx x≅ and stop.

Example 4.14 _________________________________________________________

Find the value of t where 2te t− = . This is an example of a problem where the method

of steepest descent, or another kind of optimization method can be employed, if we covert the

given problem to an optimization problem. If it exists, we can find the solution, call it optt , by

minimizing 2 2( ) ( )tf t e t−= − , where ( ) 0f t ≥ for all t , and the minimum value is ( ) 0optf t = .

The square of 2( )te t− − is used to ensure that ( )f t is nonnegative. With

2′( ) 2( )( 2 )t tf t e t e t− −= − − − , Prog. 4.8 finds optt . The results follow the program.

% Application of the method of steepest descent to solve the nonlinear

% algebraic equation: exp(-t) = t^2

clear all; clc; format long e

precision = 1e-8; % desired accuracy

mu = 1e-4; % step size

t_old = 1.0; % initial guess of the optimizer

f_of_t_initial_guess = (exp(-t_old)-t_old^2)^2

f_prime = 2*(exp(-t_old)-t_old^2)*(-exp(-t_old)-2*t_old);

t_new = t_old – mu*f_prime;

n_iterations = 1; % initialize an iterations counter

max_iter = 1e5; err = 0; % limit the number of iterations

while abs(t_new – t_old) >= precision % check for convergence

n_iterations = n_iterations + 1; % count iterations of algorithm

if n_iterations > max_iter, err = 1; break; end

t_old = t_new; % update recursion

% get new solution

f_prime = 2*(exp(-t_old)-t_old^2)*(-exp(-t_old)-2*t_old);

t_new = t_old – mu*f_prime;

end

if err == 0

t_opt = t_new

f_of_t_opt = (exp(-t_opt)-t_opt^2)^2

fprintf(‘Algorithm required: %6i iterations to converge\n’,n_iterations)

else

fprintf(‘Algorithm exceeded: %6i iterations\n’,max_iter)

end

Program 4.8. Application of the method of steepest descent.

f_of_t_initial_guess = 3.995764008937280e-001

t_opt = 7.034812353615486e-001

f_of_t_opt = 6.900847459766987e-010

Algorithm required: 13347 iterations to converge

_____________________________________________________________

The method of steepest descent has some drawbacks. It is slow to converge, which can

be improved by increasing the step size. It requires a derivative of the function to be

optimized, which for some problems may be practically impossible to obtain. However, it and

its variations are widely applied, because it is a simple algorithm to implement. The required

derivative could be approximated with, for example

( / 2) ( / 2)

‘ ( ) n nn

f x h f x h

f x

h

+

− −

(4.11)

where h is some small positive number, for example, h µ= . Then, the method requires more

function evaluations, which can increase computing time. Also, once nx is close to the

optimizer, (4.11) will give small and erratic ‘( )nf x values, making further convergence unlikely.

However, this depends on the behavior of ( )f x and the desired precision.

The MATLAB built-in function fminsearch, which is a function function, is a derivative

free optimization method [1], which makes it attractive to use. A syntax is given by

[x_opt,f_opt] = fminsearch(fun,x_initial_guess)

where fun is a function handle. For example, fun could be @myfun, where myfun is the name

of a function m-file, or it could be a handle of an anonymous function. Use doc fminsearch for

more syntax options.

Let us apply fminsearch to the problem of Example 4.14. An anonymous function is

given to define ( )f t , as follows.

>> format long e

>> f_of_t = @(t)(exp(-t)-t^2)^2; % define an anonymous function

>> f_of_t(1)

ans = 3.995764008937280e-001

>> [t_opt f_opt] = fminsearch(f_of_t,1)

t_opt = 7.035156249999998e-001

f_opt = 8.403998856470625e-009

While fminsearch is easy to use and can give good results, it can also require many

function evaluations or never converge to useful results. This depends on the behavior of the

function to be optimized and the initial guess. The function ( )f t of Example 4.14 has only one

minimum, the global minimum, and both methods of optimization used here worked very well,

regardless of the initial guess. However, there are no known methods of optimization that can

with certainty find the global minimum of a function. This depends on the behavior of the

function, and it is important to provide an initial guess in the vicinity of the optimizer.

4.7 Numerical Integration

Given is some function ( )x t that is finite for t−∞ < < + ∞ . To find ( )y t given by

( ) ( ) ( ) ( ) ,

t a t

a

y t x d x d x d t aτ τ τ τ τ τ

−∞ −∞

= = + >∫ ∫ ∫ (4.12)

is a common problem. We can write (4.12) as

( ) ( ) ( ),

t

a

y t x d y a t aτ τ= + ≥∫ (4.13)

where ( )y a is called the integration constant. If t a= , then (4.13) becomes ( ) ( )y t a y a= = .

More specifically and for example, recall Fig. 2.2, in Chapter 2, where

( )

( )

d q t

i t

d t

=

If we integrate this equation, we get

( ) ( ) ( ) ( ) ( ) ( ) ,

t t a t

a a

q t i d i d i d i d q a t aτ τ τ τ τ τ τ τ

−∞ −∞

= = + = + ≥∫ ∫ ∫ ∫ (4.14)

The second term ( )q a is the amount of charge that has passed through the wire cross-section

from t =− ∞ until t a= , and the first term is the additional amount of charge that has passed

through the wire cross-section from t a= until sometime t , where t a> .

Depending on ( )x t in (4.13), it may be difficult, if not impossible, to find a function for

( )y t , and we must resort to finding ( )y t for a set of discrete time points with some numerical

integration method. Sometimes we do not know ( )x t for all t . Instead, we only know ( )x t for

a set of discrete time points.

Without loss of generality, let 0a = , and set the upper integration limit to 0t T= .

Consider the integral

0

0

( )

T

A x t dt= ∫ (4.15)

where A is the area under the curve ( )x t from 0t = until 0t T= .

4.7.1 Euler’s Method

Let us segment the time range 0[ 0, ]T into N sub-intervals, where the width of each sub-

interval is 0 /T T N= . This is depicted in Fig. 4.9. The value of ( )x t at the beginning of each

sub-interval is ( ), 0, 1, , 1x nT n N= − . Euler’s method approximates the area under the

curve in each sub-interval by a rectangular area given by the width T times the value of ( )x t at

the beginning of the sub-interval. The integral of (4.15) is approximated by summing all of the

rectangular sub-interval areas to get

1 1

0 0

( ) ( )

N N

k k

A x kT T T x kT

− −

= =

=∑ ∑ (4.16)

where (( 1) )x N T− is the value of ( )x t at the beginning of the last sub-interval. The

approximation gets better as N is increased, which makes T smaller.

( )x t

tT 2T NT

( 1)N T−

(0)x

( )x T

(2 )x T

(( 1) )x N T−

0( )x T NT=

(4 )x T

Figure 4.9. Segmentation of ( )x t over 0[0, ]T NT= .

In view of (4.16), let

1

0

( ) ( )

n

k

S n x k T

−

=

= ∑

and then ( )A T S N . For some n , assume that we have ( ), 0, , 1x kT k n= − , which means

that we can find ( )S n . When ( )x nT becomes available, we can update ( )S n to get ( 1)S n +

with

( 1) ( ) ( ), 0, 1, , 1, (0) 0S n S n x n T n N S+ = + = − = (4.17)

This is a recursion in ( )S n , where

(1) (0) (0 ), set S(0) = 0

S(2) = S(1) + ( )

(3) (2) (2 )

( ) ( 1) (( 1) )

S S x T

x T

S S x T

S N S N x N T

= +

= +

= − + −

Each time a new value ( )x nT becomes available, ( )S n is updated to obtain ( 1)S n + . The

recursion (4.17) is operated for 0, 1, , ( 1)n N= − , with (0) 0S = .

Example 4.15 _____________________________________________________

Let ( )x t be the function defined by

1, 0

( ) sin( )

, 0

t

x t t

t

t

π

π

=

=

≠

(4.18)

This function is called the sinc function, i.e., ( ) sin c( )x t t= , and it is a built-in MATLAB function.

With

>> t = linspace(-10,10,1000); % specify time points

>> x = sinc(t); % evaluate x(t) = sin(π t)/π t at the discrete time points given by the elements of t

>> plot(t,x); grid on

>> xlabel(‘t – secs’); ylabel(‘sinc(t)’);

we get the plot shown in Fig. 4.10. The sinc function oscillates indefinitely in both directions

with an amplitude that decays to zero. The integral of ( ) sin c( )x t t= from −∞ to +∞ is one.

-10 -8 -6 -4 -2 0 2 4 6 8 10

-0.5

0

0.5

1

t – secs

si

nc

(t

)

Figure 4.10. The sinc function plotted from t=-10 to t=10.

An accurate value for the integral of the sinc function from 0t = until 10t = is:

0.4898881. Therefore, the integral from 10t = to t →+ ∞ is (0.5 – 0.4898881). Prog. 4.9 uses

Euler’s method to find this integral, and the results for N = 50, 550, 1050 follow.

% Program to find the integral of a function using Euler’s method

clear all; clc

T0 = 10; % time range is from 0 to T0

for m=1:3

N = 50+500*(m-1); % number of sub-intervals

T = T0/N; % sub-interval width

S = 0; % initialize the sum

for n = 0:N-1 % loop to include each sub-interval

t = n*T; % time point

S = S + sinc(t); % recursion

end

Area(m) = T*S; % area

N_subs(m) = N;

end

format short

table = [N_subs;Area];

disp(‘ N Area’)

fprintf(‘ %i %f \n’,table)

Program 4.9. Program to find the integral of sinc( )t using Euler’s method.

N Area

50 0.590224

550 0.498982

1050 0.494651

These results show that N must be large to achieve a small numerical integration error.

Nevertheless, depending on the behavior of ( )x t , Euler’s method is often used because of its

simplicity.

______________________________________________________________

4.7.2 Trapezoidal Rule

The trapezoidal rule is a slight modification of Euler’s method, and gives better results. Instead

of using the value of ( )x t at the beginning of each sub-interval, the trapezoidal rule uses the

average of the values of ( )x t at the beginning and end of each sub-interval, and (4.16) and

(4.17) become

1 1

0 0

( ) (( 1) )

( ) (( 1) )

2 2

( 1) ( ) ( ) (( 1) ), 0, 1, , 1, (0) 0

N N

k k

x kT x k T T

A T x kT x k T

S n S n x n T x n T n N S

− −

= =

+ +

= + +

+ = + + + = − =

∑ ∑

(4.19)

where ( / 2) ( )A T S N .

Example 4.16 _____________________________________________________

To apply the trapezoidal rule to the integration problem of Example 4.15, the inner for

loop in Prog. 4.9 is replaced by

xnT = sinc(0); % this is x(nT) for n=0 and it will be updated in the loop

for n = 0:N-1 % loop to include each sub-interval

t = (n+1)*T; % next time point

% evaluating the sinc function just once for each loop iteration

xn1T = sinc(t);

S = S + xnT + xn1T; % trapezoidal rule recursion

xnT = xn1T; % update x(nT)

end

and the area is computed with

Area(m) = (T/2)*S;

The results are:

N Area

50 0.490224

550 0.489891

1050 0.489889

Here, the result for 50N = are better than the result for 550N = using Euler’s method.

______________________________________________________________

4.7.3 Built-In Integration Functions

MATLAB has several sophisticated built-in functions for evaluating a definite integral. One of

these functions is quad. It is a function function, and the syntax to use this function is

area = quad(@my_function,a,b, tol)

where the name of the m-file that defines the function to be integrated from a to b must be

my_function.m. The fourth argument specifies an error tolerance. If tol is not included in the

input argument list, then MATLAB uses a default value of 1e-6.

The function quad was used to obtain the accurate result from integrating the sinc

function in Example 4.15. The statement that was used is: area=quad(@sinc,0,10.0,1e-8).

When MATLAB executes the quad function algorithm, it will want to pass to the user-

defined function a vector input. Therefore, the user-defined function must be properly

vectorized to return a vector. See Table 3.3 for more details about vectorizing expressions

within a function that involve function input vectors.

Example 4.17 _____________________________________________________

Let us find the amount of charge that has passed through a diode over one cycle of a

sinusoidal voltage across it. Recall the diode i-v characteristic, for which a function is given

below.

function current = diode(t)

% current through a diode when the voltage across it is v=Asin(wt)

global w A % frequency and amplitude

v = A*sin(w*t);

I_S = 1e-12; % saturation current

V_T = 25.85e-3; % thermal voltage at 300 K

current = I_S*(exp(v/V_T) – 1);

Here, the function input variable is t , which is used as the input time argument for the

sinusoidal voltage. The frequency and amplitude of the voltage are passed to this function as

global variables. This function was saved as diode.m. To see how charge moves through the

diode over one cycle of the voltage, let us segment one cycle into N sub-intervals, and

integrate the current over each sub-interval. Prog. 4.10 gives a script to plot charge movement

over one cycle.

% Program to plot charge movement through a diode over one cycle of

% a sinusoidal voltage across it

clear all; clc

global w A

frequency = 60; % specify frequency

w = 2*pi*frequency;

A = 0.8; % amplitude of sinusoid

T0 = 1/frequency; % time of one cycle

N = 100; T = T0/N; % will use N+1 plot points

tol = 1e-8; % specify tolerance

q(1) = 0; % initial charge is zero

v(1) = 0; % initial voltage value

t = 0.0:T:T0; % specify N+1 time points

for k = 1:N

v(k+1) = A*sin(w*t(k+1));

q(k+1) = q(k) + quad(@diode,t(k),t(k+1),tol); % accumulate charge

end

plotyy(t,q,t,v); grid on % use right y-axis for v

xlabel(‘time secs’); ylabel(‘Coulombs’)

title(‘Charge movement through a diode’)

Program 4.10. Script to find the amount of charge that has passed through a diode.

0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018

0

0.005

0.01

0.015

0.02

0.025

0.03

0.035

0.04

time secs

Co

ulo

m

bs

Charge movement through a diode

0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018

–

0.8

–

0.6

–

0.4

–

0.2

0

0.2

0.4

0.6

0.8

Figure 4.11. Charge movement through a diode over one voltage cycle.

Fig. 4.11 shows how charge moves through a diode over one cycle of the voltage across

it. Notice that as the voltage increases, there is no appreciable movement of charge through

the diode until the voltage reaches a threshold voltage of about 0.7 volts. Then, once the

voltage exceeds this threshold, the amount of charge that has moved through the diode

increases quickly until the voltage again goes below the threshold, after which the net amount

of charge through the diode no longer increases.

______________________________________________________________

There are several other built-in functions for numerical integration. These functions

adaptively adjust the way they work to achieve prescribed error tolerances. For example, quadl

may be more efficient than quad, while also giving more accurate results, especially for

smoothly behaving integrands. The built-in function quad2d numerically obtains double

integrals over a planar surface. See the MATLAB help facility for details and additional built-in

functions for numerical integration.

4.8 Switch – Case – Otherwise

The switch-case-otherwise structure provides a way to select which block of statements in a

set of blocks of statements should be executed. The syntax of this structure is

switch expression

case expression_1

block #1 of statements

case expression_2

block #2 of statements

insert as many additional case expression and block of statements as needed

otherwise

block #N of statements

end

where each switch key word must be matched by an end key word. The otherwise and block

#N of statements is optional. The switch expression must evaluate to a scalar or a character

string. When MATLAB encounters the switch key word, it does the following

1) evaluate the switch expression and

2) if the result of evaluating the switch expression matches the case expression_1 value, then

execute block #1 statements, and continue execution after the end statement

3) if the result of evaluating the switch expression does not match the case expression_1 value,

then skip the block #1 statements to the next case

4) if the result of evaluating the switch expression matches the case expression_2 value, then

execute block #2 statements, and continue execution after the end statement

5) if the result of evaluating the switch expression does not match the case expression_2 value,

then skip the block #2 statements to the next case, and so on

6) if the result of evaluating the switch expression does not match any case expression value,

then execute the block #N statements after the otherwise key word (if there is one) and

continue execution after the end statement.

It is useful to indent each block of statements to easily see each case block of statements and

that each switch key word is matched by an end key word. Within each block of statements

there can be any of the script flow control structures, including nested switch-case-otherwise

structures.

Example 4.18 _____________________________________________________

You are approaching an intersection with a traffic light. Different actions are required

depending on the color and behavior of the traffic light. Prog. 4.11 is a demonstration of the

application of the switch-case-otherwise structure to this situation.

% Suggested action when approaching an intersection with traffic lights

% Local conditions may require other actions

clear all; clc

disp(‘Upon approaching an intersection with working traffic lights,’)

disp(‘different actions are required.’)

behavior = input(‘Enter a traffic light behavior in single quotes: ‘);

switch behavior

case ‘red’

disp(‘Come to a complete stop; and wait for green light to proceed.’)

case ‘yellow’

disp(‘Slow down; prepare to stop.’)

case ‘green’

disp(‘Proceed with caution.’)

case ‘flashing yellow’

disp(‘Slow down; proceed with caution; be ready to stop.’)

case ‘flashing red’

disp(‘Come to a complete stop; proceed with caution.’)

otherwise

disp(‘Traffic lights may be out of service.’)

disp(‘Come to a complete stop.’)

disp(‘Local conditions may require other actions.’)

end

Program 4.11. Program to suggest action when approaching an intersection with traffic lights.

______________________________________________________________

4.9 Conclusion

Testing a condition, executing alternative parts of a script and repeatedly executing a part of a

script are essential for implementing problem solution methods. All programming languages

have means to do this. You should now know how:

• to use relational operators for formulating condition tests

• to use logical operators to make compounded statements using logical parts

• the if-elesif-else structure is used to execute alternative script parts

• to make a script more robust

• the for loop structure works to repeatedly execute a block of MATLAB statements

• to find the time required to execute a block of statements

• a histogram is obtained

• to determine the probability of an event by counting

• a pseudorandom number generator can be used to simulate random events

• a median filter works

• to use the sound device

• the while loop structure works to repeatedly execute a block of statements depending

on a condition

• Euler’s method and the trapezoidal rule are used for numerical integration

• the switch-case-otherwise structure is used to execute alternative script parts

Table 4.5 gives the MATLAB functions that were introduced in this chapter. Use

MATLAB’s help facility to learn more about these built-in functions, where you will also find out

about many other related built-in functions.

Table 4.5 Built-In MATLAB Functions Introduced in This Chapter

Function Brief Explanation

isequal(A,B) returns logical 1, meaning true, if arrays A and B are the same size

and contain the same values, and logical 0, meaning false,

otherwise

L = and(A,B) an element of the output array is set to 1 if both input arrays

contain a non-zero element at that same array location, and

otherwise, that element is set to 0

L = or(A,B) an element of the output array is set to 1 if either input array

contains a non-zero element at that same array location, and

otherwise, that element is set to 0

L = xor(A,B) an element of the output array is set to 1 if either input array

contains a non-zero element at that same array location, but not

both, and 0 if both elements are zero or non-zero

L = not(A) an element of the output array is set to 1 if A contains a zero

value

element at that same array location, and 0 otherwise

&& short logical and

|| short logical or

I = find(A) returns the linear indices corresponding to the non-zero entries of

the array A

L = sign(A) signum function, returns 1 if the element is greater than zero, 0 if

it equals zero and -1 if it is less than zero; for the non-zero

elements of complex A, sign(A) = A ./ abs(A)

tic start a stopwatch timer

T = toc saves the elapsed time in T (seconds) since the last execution of

tic

pause(t) pauses for t seconds before continuing, where t can also be a

fraction; just pause causes a script or function to stop and wait for

the user to strike any key before continuing

bar(X,Y,width) draws the columns of the M × N matrix Y as M groups of N vertical

bars; bar(Y) uses the default value of X=1:M. For vector inputs,

bar(X,Y) or bar(Y) draws length(Y) bars. There is no space

between bars if width=1.

continue pass control to the next iteration of for or while loop

break terminate execution of for or while loop

return causes a return to the invoking script or function or to the

keyboard; also terminates the keyboard mode

keyboard when placed in an m-file, stops execution of the file and gives

control to the user’s keyboard. The special status is indicated by a

K appearing before the prompt. Variables may be examined or

changed; all MATLAB commands are valid. The keyboard mode is

terminated by executing the command return.

R = rand(N,M) returns an N× M matrix containing pseudorandom numbers

uniformly distributed over the open interval (0,1)

I = randi(N,M,imax) returns an N× M matrix containing pseudorandom integer values

uniformly distributed over the range 1:imax

R = randn(N,M) returns an N× M matrix containing pseudorandom values from

the standard normal distribution

Y = mean(X) Y is a row vector containing the mean of each column of X

Y = var(X) Y is a row vector containing the variance of each column of X

I = hist(X,N) bins the elements of X into N bins

fprintf write formatted data to a display or a file; more about this later

y = median(X) returns the median value of the elements in the vector X

sound(X,fs) sends the signal in vector X, with sampling frequency fs to the

sound device; X can be a two column matrix for stereo. Clips

values outside the range -1 to +1.

soundsc(X,fs) same as sound, but autoscales

wavwrite(X,fs,’file.wav’) writes data X to a Windows WAV file specified by the file name

file.wav, with a sampling rate fs Hz; stereo data must be specified

as a matrix with two columns

Y = sinc(X) returns in the matrix Y elements found with sin(pi*x)/(pi*x) of the

elements in the matrix X

quad(func,a,b) returns an approximation of the integral from a to b of the

function specified by the function handle func

plotyy(x1,y1,x2,y2,plot_func) plots y1 versus x1 with y-axis labeling on the left and plots y2

versus x2 with y-axis labeling on the right, using the plotting

function specified by the function handle plot_func

quadl same as quad, but uses a more accurate numerical integration

method

quad2d(func,a,b,c,d) returns an approximation of the double integral of a function

f(x,y) specified by the function handle func, where ≤ ≤a x b and

≤ ≤c y d . See MATLAB help for more options.

The essentials of the MATLAB programming language have been presented in the first

four chapters of this book. Use demo to watch an excellent audio/video tutorial about if-elseif-

else, for and while loop control structures. In these chapters some special features, for

example, plotting, were introduced to make examples more complete with visual displays. In

the remaining chapters, special, unique and particularly useful features of this programming

language will be discussed.

In the next chapter the concept of a logical variable will be extended to understand its

role in Boolean algebra, logic circuits and binary arithmetic.

References

[1] Nelder, J.A. and Mead, R. (1965), “A simplex method for function minimization”, Comput. J.,

7, pp. 308–313.

[2] Li, C., Priemer, R. and Cheng, K.H. (2004), “Optimization by Random Search with Jumps”, Int.

J. for Numerical Methods in Engineering, vol. 60, pp. 1301-1315.

Problems

In all of the scripts that you are required to write for the following problems, include comments that

state the purpose and activity of the script.

Section 4.1

1) Use a = -2, b = 3 and c = 5, and find the result of each of the following relational expressions. (a)

d = a^2 >= c, (b) e = (a > b -4)+c <= b^2, (c) f = a+b == c-4.
2) For t = -1:0.1:2, w = 2*pi and x = cos(w*t), use a logical operation to find a vector y, which has
the same dimension as x, where all of the positive elements of x have been replaced by zero.
3) A = [-1 9 2; 0 -3 4], B = [], C = ones(3,2)'. Give the result of: (a) isempty(A), (b) isempty(B), (c)
logical(A), (d) any(A), (e) islogical(C), (f) C < false(2,3) | A, (g) isinteger(sqrt(A)), (h) isfloat(C).
Section 4.2
4) For a = 2, b = 0 and c = -5, give the results of: (a) d = a | (b & c), (b) e = and(a,b)|c,
(c) f = a & c | ( b & ~a), (d) g = a & c | b & ~a, (e) h = (b && c) | (b || a), (f) p = a & xor(~b,c).
5) (a) Describe the difference between (a & b) and (a && b).
(b) Describe the difference between (a | b) and (a || b).

6) For A = [0 1; 1 0], C = [0 0; 1 0] and C = A & B, give all possible B.

7) For a = -2, b = 4 and c = 2, give the results of: (a) d = (a ~= c -2) & a^2 > b,

(b) e = sqrt(a^2 + b^2) < abs(c) | (a+b)/2 >= c ,

(c) f = ~isempty(find([a b 2*c sqrt(a^2 + c^2)] > 4)).

8) For A = [-1 0.5 5; 2 -4 0; 0 0 1], (a) give y = find(A <= -1 | A >= 1) , and explain the meaning of

y. (b) let B = reshape(A,1,9) and give C = B(y).

Section 4.3

9) (a) If the inputs to Prog. 4.1 are: a = 0, b = 5, N = 9 and x = 5.6, then what is the script output?

(b) If the inputs are the same as in part (a), except that x = 2.7, then what is the script output?

10) (a) Convert Prog. 4.1 into a function, call it bin_number. The function input arguments are: a,

b, x and N, and the output arguments are: k and err. Include comments that explain the

purpose of the function and the meaning of the err codes. Provide a copy of the function.

(b) Write a script that demonstrates using the function bin_number. Provide a script listing

and demonstration of script operation.

11) Write a script that implements the flowchart given in Fig. P4.11. Use an if-else-end control

structure. Provide a script listing and a demonstration of script operation. Does the current

have to be a positive number?

Input a resistor value R

R > 0 ?

Display message

Terminating Program

NoYes

Input a current I

Calculate power

delivered to resistor

Display

power

End program

Figure P4.11. Flowchart to calculate power delivered to a resistor.

12) Write a script that inputs an integer N in the range 0 to 9. Your script must check that the

input is an integer and that it is within the required range. Otherwise, display an error message

for each input error type. Then, the script must determine if N is an odd or even integer and

display a message about the result. Provide a script listing and a demonstration of script

operation.

Section 4.4

13) Describe what is wrong with each of the following for statements.

(a) for k = 0:-1:10, (b) for n = 10:2:3, (c) for 2k = 0:5, (d) for Number = -1:-5. In each case,

what will MATLAB do when the for loop is completed with an end statement?

14) (a) Write a function script, call it inner_prod, that uses a for loop to find the inner product

A_dot_B of two vectors A and B with dimension 1xN. The function input arguments must be A

and B, and the function returns A_dot_B.

(b) Write a script to demonstrate the operation of your function where t = 0:T:T0-T, T = T0/N

secs, T0 = 02 /π ω secs, A = cos( 0ω t) and B = cos(K 0ω t) for K = 1, 2 and 5. You choose N and

0ω . For each K, what is A_dot_B? What do you think A_dot_B will be for any integer K ≠ 1?

(c) Repeat part (b) for B = sin(K 0ω t).

(d) When the inner product of two vectors is zero, it is said that the vectors are orthogonal.

Here, we find that this concept can be extended to functions, e.g., cos(M 0ω t) and cos(K 0ω t),

M ≠ K, are orthogonal functions over a period T0. What other sinusoidal functions are

orthogonal over a period?

15) Write a script that uses for loops to evaluate the power series given by (2.13) of an

exponential function 2( ) tx t e−= , for t = 0:0.1:2. Break out of the for loop that does the power

series computation when the next term in the series contributes less than 0.0001 % of the sum

of terms. For example, see Prog. 4.4. Provide a figure of a plot of your result. Include in the

figure a stem plot that uses the MATLAB built-in function exp. The figure should include a grid,

axes labels and a title.

16) Given is the circuit shown in Fig. P4.16.

(a) In terms of V, Rs and R, obtain an expression for the power P delivered to R.

(b) Write a script that receives as input the values of V and Rs, and makes a vector R with

element values that range from 0.1*Rs to 10*Rs in increments of 0.1*Rs. Use a for loop to

create a vector of power delivered to each element value of R. Then, plot the power P versus R.

Do this for V = 10 volts, Rs = 10 Ohms and Rs = 100 Ohms. For each case, from inspection of the

plots, give a value of R such that V delivers the maximum power to R. Provide a script listing

and plots, including axes labels and titles.

(c) Vectorize your expression for P, and repeat part(b) without using a for loop.

(d) In general, for a given value of Rs, what value of R enables the resistor R to receive the most

power from the voltage source V?

V R

R s

Figure P4.16. Battery with internal resistance Rs connected to a load R.

17) Write a script that starts by obtaining a vector x of N = 1e6 normally distributed random

numbers with the function randn.

(a) Obtain the mean, variance and histogram of x. Use 51 bins for the histogram. Use the

MATLAB built-in functions mean, var and hist to obtain these results.

(b) Let y = ax +b, where a = 3 and b = 2. Use a for loop to find the mean of y. Then, use another

for loop to implement (4.6). How does your result compare with using the MATLAB function var

to find the variance of y?

(c) Obtain a histogram of y, using 51 bins. How are the parameters a and b related to the

variance and mean of y?

18) Suppose that an experiment consists of tossing two die. Each die will show an integer from 1 to

6. This can be simulated with die_1 = ceil(6*rand) and die_2 = ceil(6*rand), where each

invocation of rand produces a random number from 0 to 1. Write a script that uses a for loop to

do N_exps = 1e7 experiments, and within the loop uses an if-elseif-end structure to count the

number of times die_1 + die_2 = 7 or 11. Begin your script with the statement: rng(‘default’) to

bring the random number generator to its default initial state. By counting, find the probability

that an experiment outcome will be 7 or 11.

19) Let us study the performance of the algorithm given by

1

( ) ( ( ) 2 (( 1) ) (( 2) )

4

y k T x k T x k T x k T= + − + −

for 0, 1,k = , where ( )x k T comes from sampling a sound given by

1 2( ) sin( ) 0.5 cos( )s t t tω ω= + that is corrupted by additive noise ( )n t . Let 1 2 (262)ω π=

and 2 2 (440)ω π= . The sampling rate is 8000sf = samples/second, and 1 / sT f= . The

data ( )x k T for this simulation is given by ( ) ( ) ( )x kT s kT a n kT= + , where ( )n kT is found

by using a normal random number generator and a controls the noise amplitude. To help you

get started, use

% This is an example of digital signal processing

% Program to implement a digital low-pass filter

clear all; clc;

fs =8000; T = 1/fs;

w1 = 2*pi*262; w2 = 2*pi*440;

a = 0.5; % used to see how noise affects the result

k = 0; t = k*T; time(k+1) = t; % initialize sample time

x(k+1) = sin(w1*t) + 0.5*cos(w2*t) + a*randn; % first sample

% since there is no x(-T) and x(-2T) use

y(k+1) = x(k+1)/4;

k = k+1; t = k*T; time(k+1) = t; % save time for plotting

x(k+1) = sin(w1*t) + 0.5*cos(w2*t) + a*randn; % second sample

% since there is no x(-T) use

y(k+1) = (x(k+1)+2*x(k))/4;

N = 100; % process N-2 more samples of the input

for k = 2:N-1

t = k*T; time(k+1) = t; % save time for plotting

x(k+1) = sin(w1*t) + 0.5*cos(w2*t) + a*randn; % kth sample

y(k+1) =

Complete the for loop with the given algorithm, and then plot x and y versus time. What does

the algorithm appear to do?

Section 4.5

20) Do Prob. 4.11 using a while loop. The script should allow for repeated execution of the script

until a program user enters an R value less than or equal to zero. For example,

clear all; clc;

while 1

R = input(‘Enter a resistance value: ‘)

if R <= 0
disp('Terminating program')
break
end
Complete the script. For a nice looking output, use, for example, fprintf('Power delivered =
%6.2f watts\n', P). Provide a script listing, and demonstrate its operation.
21) Do Prob. 4.15 using a while loop inside a for loop. With the for loop index ranging from n
= 1 to n = N, where N = length(t), then, for example, the while loop could be
x = -2.0*t(n);
exp_sum = 1.0; % first term in power series
next_term = x; k = 1; % second term in power series
while abs(next_term) > 1e-6*abs(exp_sum)

exp_sum = exp_sum + next_term;

k = k+1;

next_term = next_term*x/k;

end

exp_func(n) = exp_sum + next_term; % use last term

Explain why the abs function must be used. Complete the script, and provide results as

described in Prob. 4.15.

Section 4.6

22) Consider again Prob. 4.16. Write a script that applies the method of steepest ascent, (4.10), to

find the value of R that maximizes the power delivered to it by the voltage source V. Use input

statements to obtain the values of V and Rs, and a while loop to implement an update method

to find Rnew given Rold similar to the update method of (4.9). See Example 4.14. Find a suitable

value of the step size by trial and error. Use Rold = Rs/2 as the initial guess of Ropt. Provide a

script listing and a demonstration of its operation.

23) Given a number x, find y = x1/2. Of course, y = sqrt(x). However, an algorithm must be

designed to find y. This can be done by converting the given problem to an optimization

problem. Consider minimizing the function f(y) = (y2 – x)2, where f ‘(y) = 2(y2-x)(2y). Write a

script that uses a while loop within a while loop to apply the method of steepest descent. The

outer while loop must contain an input statement to receive a value of x, and if the input is not

greater than zero, then break out of this loop, which terminates the script. If the input x is

greater than zero, then the inner while loop must apply the method of steepest descent to find

y. Implement an update method to find ynew given yold like the update method of (4.9). Use yold

= x as an initial guess of yopt. See Example 4.14. Find a suitable step size by trial and error.

Provide a script listing and a demonstration of its operation for three positive values of x.

24) Repeat Prob. 4.23, but replace the inner while loop with an invocation of the MATLAB built-in

function fminsearch. Create a function m-file to define the function to be optimized by

fminsearch. Provide listings of your script, the function m-file, and a demonstration of script

operation.

Section 4.7

25) Given is the function x(t) = e-at – e-bt. Let a = 1 and b = 2.

(a) Manually, find an expression for the area A under x(t) from t = 0 to t = 5.

Write a script to:

(b) plot x(t) for t = 0 to t = 5. Use N = 501 values of t in the given time range

(c) calculate A using the expression found in part(a)

(d) apply Euler’s method to find an approximation of A for N = 11, 111 and 211

(e) output a table of the difference between A and its approximation for each N.

Use format long e for all script output. Provide a listing of your script and the results.

26) Repeat Prob. 4.25, but use the trapezoidal rule.

27) For the x(t) given in Prob. 4.25, use the built-in MATLAB function quad to find the area A

under x(t) from t = 0 to t = 5. Use an anonymous function to define x(t). Compare the result to

A found in Prob. 4.25, part(a). Use format long e for all script output. Provide a listing of

your script and results.

Section 4.8

28) Write a function script, call it linear_eqs_solver, where the inputs are a square matrix A, a

column vector Y and a character string named solution_method, and the outputs are the

condition number K of A and the solution X of AX=Y, assuming that the determinant of A is not

zero. First, check the determinant of A, and if it is zero, set solution_method = ‘none’. The

possible solution methods you should use are: inv (the function), left matrix divide, eigenvalues

and eigenvectors or singular value decomposition, as designated with solution_method = ‘inv’,

‘left_mat_div’, ‘eigen’ or ‘singular_decomp’. Use the switch-case-otherwise-end structure to

select the method to be used. Use the otherwise option to return an empty vector X and K = inf

in the event that the character string solution_method is not one of the allowed solution

method options. Write a script that demonstrates the application of your function. Provide

script and function listings and script output.

29) Repeat Prob. 4.18, where within the for loop a switch-case-end structure is used to count the

number of times that an experiment outcome is 7 or 11. Use: die_1 + die_2, for the switch

expression and case 7 and case 11.

- Table 4.5 gives the MATLAB functions that were introduced in this chapter. Use MATLAB’s help facility to learn more about these built-in functions, where you will also find out about many other related built-in functions.
- The essentials of the MATLAB programming language have been presented in the first four chapters of this book. Use demo to watch an excellent audio/video tutorial about if-elseif-else, for and while loop control structures. In these chapters some s…
- In the next chapter the concept of a logical variable will be extended to understand its role in Boolean algebra, logic circuits and binary arithmetic.