# Write a Matlab Program

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

___________________________________________________________________

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
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

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

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

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
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
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

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

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.